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Preface  to  the  final  report.  -  H.  L.  Kuo,  the  principal  investagator 

As  stated  in  our  renewal  proposal,  our  tasks  during  the  first  year  of 
this  contract  are  synthesis  of  shallow  and  deep  convection  systems  and  for- 
,  mutation  of  tractable  dynamic  models  of  the  squall-line  type  thunderstorms, 

and  the  tasks  during  the  second  year  are  to  analyze  the  influences  of  the 
vertical  shear  of  the  mean  wind  on  the  organization  of  these  meso-scale 
thunderstorms  and  the  integration  of  the  squall-line  models  formulated,  plus 
investigations  of  other  subjects  related  to  severe  storms.  We  have  by  and 
large  followed  this  research  scheme  and  the  results  we  obtained  are  presented 
in  the  four  papers  included  in  this  final  report. 

It  is  well  known  that  all  squall-line  type  thunderstorms  are  prominently 
two-dimensional  in  structure  even  though  new  cells  often  form  at  the  leading 
edge  and  old  cells  decay  at  the  trailing  edge  along  the  line,  so  that  the 
variation  in  the  direction  parallel  to  the  squall-line  is  also  of  some  impor¬ 
tance.  Further,  the  velocity  in  the  direction  parallel  to  the  squall-line 
is  often  also  of  importance  both  for  the  generation  and  for  the  development 
of  the  storm  system  and  it  usually  varies  significantly  both  in  the  direction 
parallel  to  the  line  and  perpendicular  to  the  line,  as  examplified  by  the 
moisture  leaden  southerly  current  from  the  gulf  of  Mexico  for  most  of  the 
squall-line  type  storms  in  mid-west  United  States.  In  order  to  include 
these  effects  and  yet  avoid  the  use  of  the  very  complicated  and  time  consuming 

*  A 

three-dimensional  model,  a  multi-two-dimensional  model  has  been  developed  by 
»  the  principal  investigator  in  which  the  variable  x(x,y,z,t)  approximated 

by  the  sum  of  Its  projections  X)(x,z,t)  and  x2(y>z,t)  on  the  xz-  and  the  yz- 
pianes.  A  third  component  x^(x<y«t)  can  also  be  added  to  the  system  to 


A. 


represent  the  prominent  z  -  independent  part  of  the  flow  field  such  as  the 
rotational  motion  in  tornadoes.  We  shall  use  this  model  to  investigate  the 
various  aspects  of  the  squall-line  type  storms,  especially  those  aspects  • 

which  are  related  to  the  development  and  the  structure  of  these  storms  and 

a 

the  formation  of  new  disturbances.  This  model  and  the  results  obtained  from 
it  up  to  now  are  presented  in  paper  (1). 

One  very  characteristic  and  dynamically  also  very  important  feature  of 
the  squall-line  type  thunderstorm  is  that  the  updraft  and  downdraft  usually 
slope  upward  in  the  direction  opposite  to  the  vertical  shear  of  the  mean 
wind.  Even  though  this  aspect  of  the  squall-line  structure  has  been 
discussed  by  many  dynamic  meteorologist,  no  satisfactory  explanation  has 
yet  been  found  for  it.  Here  we  shall  consider  it  as  one  of  our  major  goals 
of  this  research.  Our  analyses  on  the  various  factors  which  influence  the 
vertical  structure  of  these  prominently  two-dimensional  squall-line  type 
storms  indicate  that  this  arrangement  is  the  result  of  water  loading  and 
precipitation  because,  with  the  upshear  slope,  the  precipitation  process  relieves  the 
liquid  and  solid  water  loading  from  the  updraft  and  turns  the  downward  drag 
exerted  by  them  to  an  accelerating  force  for  the  downdraft,  thereby  making 
the  storm  system  most  efficient  energetically  whereas  with  either  a  vertical 
or  a  downshear  slope  the  downward  drag  always  tends  to  destroy  the  updraft. 

It  is  also  evident  that  a  sloped  updraft-downdraft  can  be  maintained  in  convec¬ 
tive  systems  only  when  a  vertically  shearing  mean  wind  is  present.  This  r 

conclusion  can  also  be  reached  from  their  influences  on  the  vorticity  field. 

We  consider  the  verification  of  this  hypothesis  as  one  of  our  major  goals  of 
our  research,  and  shall  treat  it  both  as  an  integral  part  of  the  general  squall¬ 
line  development  problem  by  the  use  of  our  mul ti -two-dimensional  squall -line 


model  in  work  (1)  in  the  future  and  directly  as  a  separate  problem  by  the 
use  of  the  much  simpler  purely  two-dimensional  model.  This  latter  approach 
is  more  advantages  so  far  as  the  upshear  slope  problem  is  concerned  because 
it  addresses  only  to  this  problem  and  we  have  adopted  it  in  paper  (2).  In 
this  approach,  we  start  from  a  finite  amplitude  disturbance  with  updraft  and 
downdraft  adjacent  to  each  other  and  let  It  evolve  into  a  steady  state  with 
the  influences  of  liquid  and  solid  water  loading  and  evaporation  included, 
and  try  to  find  out  whether  an  upshear  slope  will  result.  A  case  study  of 
a  real  squall-line  type  storm  is  also  included  in  this  research  both  for 
demonstrating  the  nearly  two-dimensional  structure  of  the  real  storm  and  to 
illustrate  the  unstable  stratification  of  the  atmosphere  under  which  the 
storm  develops. 

As  in  all  meso-scale  disturbance  modellings,  we  always  have  to  deal 
with  the  boundary  conditions  at  the  side  boundaries.  Since  no  real  boundary 
is  present, any  boundary  condition  imposed  may  generate  some  ficticious 
responses.  It  is  thought  that  the  radiative  boundary  condition  which  allows 
outward  propagating  disturbances  to  go  out  un-impedded  will  do  least  damage 
to  the  results.  Paper  (3)  is  concerned  with  the  most  advantageous  treatment 
of  the  conditions  at  these  open  boundaries. 

Aside  from  these  works  which  deal  either  with  the  dynamics  of  the  squall - 
line  type  disturbance  directly  or  is  concerned  with  the  most  advantageous 
treatment  of  the  open  boundary  conditions  used  in  the  modeling,  Dr.  Raymond 
and  I  have  also  worked  on  higher  order  similarity  solution  for  the  tornado¬ 
like  vortex,  which  represents  an  extension  of  my  first  order  similarity 
solution  and  allows  the  inclusion  of  the  influence  of  variation  of  the 
horizontal  velocities  with  height  on  the  vortex  structure.  The  results  of 
this  study  is  presented  In  paper  (4). 
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Multl-Two-DImensional  Model  of  Squall  Line  Type  Disturbances 


Abstract 

A  multi -two  dimensional  dynamic  model  has  been  formulated  for  the 
simulation  of  the  development  of  squall  line  type  thunderstorms  by 
approximating  every  flow  variable  by  the  sum  of  three  components  defined 
by  averaging  in  y,  x  and  z  directions  respectively,  so  that  every  compo¬ 
nent  taken  by  itself  represents  a  purely  two-dimensional  system.  This 
system  of  equations  is  integrated  from  an  initial  state  characterized  by 
the  geostrophic  basic  currents  U(y,z)  and  V(x,z)  to  24  hours  under  four 
different  conditions,  namely,  cases  la  and  2a  are  from  the  dry  and  wet 
versions  of  the  purely  two-dimensional  model  while  cases  lb  and  2b  are 
from  the  dry  and  wet  versions  of  the  coupled  two-dimensional  model.  It 
is  found  that  in  general  the  disturbance  evolves  to  a  cold  front  type, 
and  the  variation  of  the  flow  variables  in  the  direction  parallel  to  the 
squall  line  enhances  the  intensity.  Further,  the  vertical  velocities  in 
the  xz-plane  usually  show  a  double  maxima,  one  at  the  2.5  km  level  and 
another  at  7  km,  while  in  the  yz-plane  the  disturbance  consists  of  a 
number  of  cells  with  a  horizontal  scale  of  about  400  km. 
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1.  Introduction 

It  is  well  known  that  the  structure  of  the  squall  line  is  predominantly 
two-dimensional  but  some  of  the  very  important  behaviors  of  the  system,  such 
as  the  formation  of  new  convection  cells  at  the  leading  edge  and  the  supply  of 
moist  air  to  the  storm,  are  also  significantly  influenced  by  the  variation  of 
the  flow  variables  in  the  direction  parallel  to  the  squall  line.  In  order  to 
take  these  effects  into  account  and  yet  avoid  the  use  of  the  very  time  consuming 
completely  three-dimensional  model,  a  multi-two  dimensional  dynamic  model  has 
been  formulated  on  the  basis  of  the  mean  flow  condition  defined  by  the  averages 
of  the  flow  variables  in  y,  x  and  z-  directions,  which  are  taken  as  parrallel 
and  perpendicular  to  the  squall  line  and  in  vertical  direction,  respectively. 

That  is  to  say,  we  replace  or  approximate  the  flow  variable  x  by  the  variable 
Xg  defined  by 

X$(x,y,z,t)  -  Xj  (x,z,t)  +  X2(y.*»t)  +  XjU.y.t),  (1) 

and  take  Xj »  X2  and  Xj  as  the  averages  of  x(*,y,z,t)  in  y-  and  x-  and  z-di rections, 
respectively,  viz., 

I  ,  -  1  2l  +  °2  1 

X  dy,  X2  -  X  -  /  X  dx,  x3  ■  X  -  h 

2  X1 

Through  this  splitting  procedure,  every  dynamic  equation  is  replaced  by  three 
seperate  equations  containing  the  two-dimensional  variables  in  the  particular 
direction  plus  average  values  of  other  quantities. 


*1 
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n 
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X  dz 


(la-c) 
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2.  The  Governing  Equations 

On  averaging  the  equations  of  motion,  the  heat  equation,  the  water  vapor 
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j  » 


and  cloud  liquid  water  mixing  ratio  equations  and  the  anelastic  continuity 
equation  in  x-,  y-  and  z-  directions  we  then  obtain  the  following  set  of 
equations 


J3 

* 
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■  o  , 

(8a) 

9PoV2  .  3pow2 

3y  3z 

1  0  , 

(8b) 

3u.  3v. 

■S-i  +  7: — —  rn  0 

3x  3y 

t 

(8c) 

where  the  diffusion  terms  involving  the  repeated  Index  1  imply  summation  over 
i  *  1,2,3.  Here  it  «  p/pQ,  *  ■  R/cp,  pQ  -  P0(z)  is  the  undisturbed  density  of 
the  air,  8  is  potential  temperature  and  9q(z)  is  (ts  undisturbed  value,  6y  ■ 

6  +  0.6lq(|-)e  is  the  virtual  potential  temperature,  q  and  c  are  the  water  vapor 
and  cloud  liquid  water  mixing  ratios,  V|  and  v#j  are  the  eddy  viscosity  and  eddy 
conduction  coefficients  in  x{  -direction,  L  is  the  latent  heat  of  condensation 
and  sublimation,  6  is  the  condensation  rate,  j  ■  1,2,3  end 
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°1*1  “  *lt  +  *U1  +  U2  +  fty  *U+  S  +  w2*  *  ]z  +  vl'  (*2  +  Xl\  ■ 


°2X2  “*2t  +  (V,>  +  ^  +  *o„+  (w<>  +  *i)  +  (X,  +X,)„  . 


'1  3  2y 


r  2z  “2  W1  3  x 


°3*3  -  X3t  +  (u3  +  0,  +  u2)  X  3x  +  (v3  ♦  v,  ♦  02)  X3y  , 


where  X  stands  for  any  one  of  the  dependent  variables  u,v,w,6  and  q  and  the 
subscripts  t,x,y,z  denote  partial  differentiations,  and  vj  ■  v^  ,  u£  ■ 
u2  "  u2*  ^or  simplicity,  we  shall  limit  ourselves  in  this  investigation  to 
the  xz-plane  and  yz-plane  systems  only  by  assuming  X *  0.  According  to  the 
continuity  equations  (8a)  and  (  8b),  the  velocities  Uj  and  w^  and  v2  and  w2 


can  be  expressed  in  terms  of  the  momentum  stream  functions  and  ^2, 


respectively,  viz.. 


U1  "  -  “o  *12*  W1  "  “O  *lx 


-  ao  *2z’  w2  -  ao  hy  ’ 


(10a) 

(10b) 


where  aQ  ■  1 /pQ -  For  convenience,  we  shall  eliminate  the  pressure  gradient 


terms  in  (2)  -  (4)  by  cross-differentiation  between  (2)  and  (4)  and  between 
(3)  and  (4)  to  obtain  the  equations  for  the  vorti cities  and  n2  in  y-  and 

x-di rections,  and  write  the  two  sets  of  prognostic  equations  as 
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where  6^  and  62  are  the  contributions  of  the  two  two-dimensional  flow  fields 

to  the  total  condensation  rate  a  =  0$  =  +  02 ,  L  i s  the  latent  heat  of  conden- 

2 

sation  and  sublimation  and  r»|  >  r^,  °2  and  the  operators,^  2  and&z  are 

given  by 
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Here  the  factor  k  in^  and^  is  introduced  to  represent  the  difference  in 
the  values  of  the  eddy  diffusion  and  eddy  conduction  coefficients  in  horizontal 
and  vertical  directions. 

The  condensation  rates  6^  and  'n  (13)  “  (20)  are  calculated  by  the 
scheme  developed  by  Kuo  and  Qian(  1981)  which  is  to  compute  the  temperature 
and  humidity  T1  and  q'  at  time  t  rAt  from  the  above  mentioned  equations  with 
6j  and  62  deleted,  and  then  determine  the  amount  of  condensation  during  the  time 
from  t  to  t  *At  when  q^,  ■  q  j  (ttAtJ+q^fAt)  is  higher  than  k  times  the  saturation 
mixing  ratio  9s(Tg)  where  *  Tj  (t*At)+T^ttAt)  and  k  is  taken  as  slightly  less 
than  1  to  take  into  account  the  situation  where  condensation  takes  place  only 
in  a  fraction  of  the  area  represented  by  the  discrete  grid  point.  Our  formula 
used  for  calculating  6j  is 
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provided  the  first  factor  is  positive.  When  this  factor  is  negative,  6^  is 

set  to  zero.  When  condensation  takes  place,  the  final  temperature  T.ft+At)  is  given *by 

J 

T'-(trAt)  +  L<5.At/C  and  the  mixing  ratio  is  given  by  q  (T - ) ,  while  the  liquid 
J  J  P  "  J  » 

water  mixing  ratio  C  is  calculated  by  (15)  and (  20).  The  detailed  scheme  of 
C-calculation  is  given  in  Raymond  and  Kuo  (t98l). 

3-  The  Initial  State 

In  this  preliminary  study  we  have  taken  the  initial  perturbation  vorticity 
fields  as  zero  and  assigned  the  following  horizontal  velocity  fields: 


Uj(x,z)  *  2  +  3*0  tanh  (z/z^) 

(23a) 

vj  (x,z)  -  -  j-  vmj  1  -  tanh [8 (x  -  az  -  xq)]  J 

+  vm  exp  j  -R-J  [(z  -  zo)2  +  Yj  (x  -  xQ)2]  j  ; 

(23b) 

u2(y,z)  -  um  exp  j-R2(Y,(z  -  z2)2  +  Y2(y  -  yQ) 23 }  , 

(24a) 

v2(y,z)  -  umexP  J-r1z  -  Zq)2}, 

(24b) 

where 

6  -  (48m)"1,  a  -  -100  ,  R,  -  2,000//2m  ,  R2  -  10,000/^m  , 

R^  *  2000/^2m  ,  zq  -  4000m  ,  Zj  -  5000m  ,  z2  ■  9,000  m  , 

xQ  -  480  km  ,  yQ  =  1 440  km  ,  Yj  -  0.03  ,  Y2  ■  0.02  ,  Yj  ■  8.0  , 

v  ■  15  m  s  '  ,  u  *25  ms*, 
m  m 

2  -1 

The  eddy  diffusion  coefficients  are  currently  taken  as  v  ■  5  m  s  , 

© 

2  “  1 

v  *  3*5  m  s  and  the  factor  K  Is  taken  as  1000  for  the  dry  case  and  400 

for  the  wet  case.  The  density  stratification  factor  is  taken  as  equal  to 
-4  - 1 

0.95x10  m  .  Here  only  a  small  vertical  shear  is  assigned  to  Uj  because  we 
wanted  to  compare  the  results  obtained  from  the  purely  two-dimensional  system 


A 
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with  X2  *  0  with  the  results  obtained  by  Olanski  and  Ross  (1977)  which  were 
based  on  a  small  shear  In  Uj  but  a  relatively  large  vertical  shear  has  been 
included  In  u2  In  the  second  system.  This  apparent  inconsistency  will  be 
removed  in  our  further  investigations. 

The  initial  temperature  fields  O^x.z)  and  6^  y,z)  are  taken  as  in  thermal 
wind  balance  with  the  initial  velocity  fields  vj  and  u2  and  hence  they  are 
obtained  from  the  following  relations: 

1_301  ,  f  9V1 

eS'ST*  g  sr  •  (25a) 

1  3e2  f  3u2 

•  (25b) 

The  initial  humidity  is  expressed  in  terms  of  the  relative  humidity  Y 
which  is  assumed  to  increase  linearly  from  its  surface  value  of  0.8  to  its 
maximum  value  of  0.90  at  the  2  km  level  and  then  decreases  linearly  with 
increasing  height  above  this  level  according  to  the  following  formula: 

y( 4  -  0.8  +  0.05  z  for  0  <  z  <  2.0  km, 

-  0.983333  -  0.04166  z  for  2  km  <  z  ,  (26) 

where  z  is  in  km.  The  humidity  field  is  taken  as  unaltered  during  the  first 
6  hours  and  the  mixing  ratio  equations  (14),  (19)  and  (19)  and  (20)  are  acti¬ 
vated  at  t  ■  6  hours. 

4.  The  Boundary  Conditions 

4a.  The  bottom  boundary  conditions  -  The  conditions  we  impose  at  the  surface 


level  z 

■  0  are 

'i'j  -  0, 

ulz-0,  e|z-0,  q|z-o; 

(27a-d) 

\P2  -  0, 

V2z-°'  e2z-0,  ,2z-0. 

(28a-d) 

The  first  two  conditions  of  these  sets  are  equivalent  to  the  vanishing  of  the 
vorticities  rj|  and  n2  at  the  surface.  The  conditions  for  v  j  and  u2  are  that 
they  satisfy  the  thermal  wind  relations  (  25a)and  (25b/  at  the  surface. 


4b.  Boundary  Conditions  At  the  Top  z  -  H 


At  the  top  we  assume  that  both  rtj2.  0j2,  V^z  and  u2z  retain  their 
initial  values  while  qj(H)  is  obtained  from  linear  extrapolation  from  below, 
and  6jz  is  taken  as  zero. 

4c.  The  Side  Boundary  Conditions 

The  side  boundaries  x  *  0  and  x  ■  I) ,  y  ■  o  and  y  »  D,  of  the  system  are 
taken  as  open  boundaries  and  the  conditions  we  impose  on  the  flow  variables  at 


these  boundaries  are  the  following: 


8nl  +  f  8nl  .  ,  3nl 

w  +  cn  +  ci2 


(29a) 


3X, 

TGT"  0  *  X1  “  vr  er  V  ci; 


(29b) 


3n, 


an. 


W~  +  c2i  3T  +  C22  3z~ 


1 - 


(30a) 


3X2 

=5T 


o  ,  X,  ■  6,,  q,,  c 


2  • 


(30b) 


Here  C^  and  C^  are  taken  as  the  x-component  and  the  z-component  of  the  phase 
velocity  of  the  vorticity  field  while  C2j  and  C22  are  the  y-component  and 
the  z-component  of  the  phase  velocity  ?2  of  n2>  respectively.  Further,  we  assume 
that  and  ?2  are  proportional  to  Vn^  and  Vn2,  respectively,  so  that  we  have 


'll 


nlt  nlx 
“2 - ~ 

nlx  +  nlz 


r  -  n2t  n2y 

21  n2  +  n2 
n2y  n2z 


C12"- 


C22  “  " 


nlt  niz 
nix  +  niz 


n2t  n2z 
n2y  +  n2z 


(31a, b) 


(32a, b)r 


The  values  of  these  phase  velocity  components  at  the  boundaries  are  obtained 

from  extrapolations  from  the  interior  and  the  generalized  multi-dimensional  radiation 

conditions  (29a)  and  (30a)  are  applied  only  when  Cj j  and  C2j  are  in  the  outward 
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directions,  whereas  they  are  set  to  zero  when  the  phase  velocities  are  inward. 

The  advantages  of  this  mul ti -dimensional  radiation  condition  and  the  ways  of 
its  implementation  have  been  demonstrated  and  described  by  the  authors  in 
another  paper  (Raymond  and  Kuo,  1981). 

5.  Scheme  of  Integration 

The  prognostic  equations  (11)  -  (20)  are  integrated  at  discrete  grid  points 
and  time  defined  by  Xj  -  iAx,  Zj  =  jAz,  yk  -  kAy,  t  =  nAt,  with  i ,  j ,  k  ranging 
from  1  to  I,  J,  K,  respectively.  Ax  ■  20  km,  Az  *  500m,  Ay  -  40  km,  and  with  I 
and  K  both  taken  as  equal  to  73  and  K  -  29  in  the  first  experiment  and  I  ■  K 
=  99,  J  ■  37  in  the  second  experiment  in  order  to  reduce  the  errors  created  by 
the  side  and  top  boundaries.  In  our  approach,  the  space  variations  are  expressed 
by  finite  element  representation  based  on  bilinear  basis  functions  while  the  time 
integration  is  by  centered  difference  leap-frog  method,  with  At  taken  as  100  sec. 

The  stream  functions  ^  and  ^  are  obtained  by  solving  the  poisson  equa¬ 
tions  (21a)  and  (21b)  together  with  the  boundary  conditions  for  these  variables. 
However,  since  the  mesoscale  disturbances  under  consideration  are  quasi-hydro¬ 
static  in  nature  and  since  our  use  of  the  large  horizontal  grid  sizes  Ax»Az 
and  Ay»Az  is  based  on  this  property,  we  can  take  the  non-hydrostatic  terms 
and  \p2yy  of  (21a)  and  (21b)  as  small  corrections  and  represent  and 

by  series  expansions  in  powers  of  the  small  parameter  e  *  ( Az/^J  ,  where 
Az/Ax  can  be  taken  as  representative  of  the  ratio  between  the  vertical  and  the 
horizontal  scales  of  variation  of  the  disturbance.  This  scheme  is  equivalent 
to  the  one  adopted  by  Orlanskl  (1981),  which  is  to  represent  by  the  following 

series 

J"  (1) 

“  J-  (33) 

*  j-o  J 

with  the  different  components  satisfying  the  following  equations 
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(1)  (1) 
Ui  +  o  \ b 

-  ponl  ’ 

ozz 

z  oz 

(1) 

(1) 

(1) 

ill,  +  a 

yJzz  z  ~j z 

"  Vlxx, 

j>l. 


(34a) 

(34b) 


0) 


Here  i|>Q  represents  the  hydrostatic  component  of  while  the  \J>j,s  ,  J>1 , 

are  the  non-hydrostatic  contributions.  In  these  forms,  the  various  components 

(1) 

are  solved  with  respect  to  z  only.  The  associated  boundary  conditions  are 

(1)  x 

4>.  -  o  at  z  ■  o  for  all  j,  I35aj 

io  o) 

if>o  ■  ¥  «  initial  value  of  i|>j,  ipj  ■  0,  j  >  I  at  z  ■  H  (35b) 

Since  e  is  rather  small  for  the  disturbances  under  consideration,  the  series 

(33)  converges  rapidly  and  taking  nj ■  2  is  sufficient  for  our  purposes. 

At  the  lateral  boundaries,  say  at  x  ■  o  and  x  **  Ik  the  right  hand  side 

0)  . 

of  (34b)  must  be  approximated,  or  with  eva  1  ua  ted  from  a  given  relation  such 

(i)  J 

as  w^  ■  o.  This  procedure  is  utilized  in  the  ^  and  ^  calculations. 

6.  Some  Preliminary  Results 

The  system  of  equations  given  above  have  been  integrated  under 
four  different  conditions  in  order  to  investigate  the  influences 


of  the  variation  of  the  flow  variables  in  y-direction  on  the  squall  line  type 
disturbance,  namely,  cases  la  and  2a  are  from  the  dry  and  wet  versions  of  the 
purely  two-dimensional  (x£  “  0)  model  while  cases  lb  and  2b  are  from  the  dry  and 
wet  versions  of  the  coupled  multi  two-dimensional  model.  Here  we  shall  present 
some  of  these  preliminary  results  and  discuss  them  briefly. 

Case  la  and  lb.  Each  of  these  two  dry  cases  has  been  integrated  to  24  hours  from 
the  initial  conditions  given  by  (23a)  -  (24b)  and  the  thermal  wind  relations 
(25a, b)  and  the  vertical  velocity  Wj(x,z)  given  by  these  two  different  versions 
at  t  ■  15  hr.  are  represented  in  Figures  1  and  2,  respectively.  From  Fig.  1  we 
see  that  the  2  units  (  ■  0.15  cm  s  *)  upward  velocity  given  by  the  purely  two- 
dimensional  xz-plane  model  reaches  to  just  the  z  ■  5*3  km  level  and  is  centered 
near  the  x  ■  700  km  mark  on  the  horizontal  scale,  and  the  maximum  upward  velocity 


t 
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In  the  xz-plane  given  by  the  coupled  mul tl -two-dimensional  model  at  t  -  15  hr. 
•reaches  the  z  *  8.67  km  level,  and  this  has  two  maxima  Instead  of  one, 
namely,  one  at  the  2.5  km  level  with  wmaxj  “  0.49  cm  s  '  and  another  one  at  the 
7.0  km  level,  with  “  0.38  cm  s  These  results  show  that  the  variations 

of  the  flow  variables  in  y-direction  contributes  significantly  to  the  flow  In 
the  xz-plane  even  under  dry  conditions  in  addition  to  their  influences  on  the 
flows  in  the  yz-plane,  and  we  expect  that  this  influence  will  become  larger 
when  condensation  process  is  included. 

Case  2a.  This  moist  purely  two-dimensional  case  has  been  integrated  for  45 
hours  and  the  results  show  that  the  vertical  motion  w^  develops  rapidly  just 
after  t  ■  6  hr.  and  the  maximum  upward  velocity  is  reached  at  t  -  16  hr.,  while 
the  intensity  of  the  disturbance  decreases  gradually  afterward.  The  distribu¬ 
tions  of  the  vertical  velocity  Wj  and  the  northward  velocity  Vj  given  by  this 
model  at  t  ■  15  hr.  in  the  xz-plane  are  illustrated  in  figures  3  and  4,  respec¬ 
tively.  From  fig.  3  we  see  that  in  this  moist  case  with  condensation  the 
maximum  vertical  velocity  reached  at  t  ■  15  hr  is  9.25  cm  s  1  which  is  about 
19.2  times  that  of  the  dry  case  la  and  it  also  occurs  at  the  3  km  level.  From 
fig.  4  we  see  that  the  positive  and  negative  northward  velocities  are  separated 
by  a  surface  which  slopes  upward  toward  west  and  therefore  the  disturbance  is 
of  the  cold  front  type,  and  the  northward  velocity  reaches  22.1  m  s  1  at  the 
5  km  level.  Here  the  vertical  circulation  in  the  xz-plane  dissipates  rapidly 
after  18  hours  but  the  yj  field  dissipates  more  slowly  as  it  is  nearly  in 
thermal  wind  balance  with  the  temperature  distribution  9j(x,z). 

Case  2b.  This  moist  case  has  been  integrated  from  the  coupled  multi -two-dimen¬ 
sional  model  up  to  24  hours  and  the  results  show  that  for  this  case  the  vertical 
motion  Wj  in  the  xz-plane  did  not  start  to  increase  until  t  ■  12  hr.  but  it 
continued  to  increase  thereafter  for  a  much  longer  time  until  t  ■  27  hr.,  and 
its  maximum  value  at  t  ■  24  hr.  is  27.5  cm  s  *  which  is  about  3  times 
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the  9*25  cm  s  *  value  In  fig. 4.  The  maximum  value  of  v-j  given  by  this 
model  at  t  ■  24  hr.  is  25.2  m  s  '  which  is  only  slightly  higher  than  the  22.1  ms 
value  given  by  2a.  The  distributions  of  and  Vj  in  the  xz-plane  for  this 
case  at  t  ■  15  hr.  and  t  ■  24  hr.  are  illustrated  in  Figs.  5  and  6  and  Figs. 

7  and  8,  respectively.  From  fig.  5  we  also  see  that  Wj  has  two  maxima  in 
the  vertical,  one  at  the  2.5  km  level  and  the  other  at  the  7  km  level.  On 
comparing  fig.  6  with  fig.  8  we  see  that  the  cold  front  has  moved  about  100  km 
in  9  hours  and  the  southerly  jet  has  also  widened  somewhat  during  this  time. 

In  addition,  we  also  present  the  stream  function  ^(y.z)  for  t  ■  24  hr.  in 
fig.  9  to  illustrate  the  pattern  of  the  meridional  circulation.  This  figure 
shows  that  this  flow  is  composed  of  a  number  of  cells  on  the  forward  northern 
edge,  with  a  horizontal  scale  of  about  400  km,  and  the  maximum  value  of  the 
vertical  velocity  w 2  of  this  circulation  is  9-32  cm  s  ' ,  which  is  about  1/3 
of  wlmax  in  fig.  8.  These  results  show  that  the  variations  of  the  variables 
in  y-direction  is  very  important  for  the  development  of  the  vertical  circula¬ 
tions  in  both  the  xz-plane  and  in  the  yz-plane  for  the  squall-line  type 
disturbance.  The  cellular  structure  of  the  circulation  in  the  yz-plane  is 
apparently  the  result  of  the  instability  of  the  vertical  shear  of  the  u2  profile 
with  regard  to  overturning  or  symmetric  perturbation  in  the  yz-plane,  and  the 
location  of  these  cells  indicates  that  new  cells  are  being  generated  on  the 
forward  edge  of  the  squat  1-1 ine. 

References 

Kuo,  H.  L.  and  Y.  F.  Qian,  1981:  Influence  of  Tibetian  plateau  on  Cumulative 
and  Diurnal  Changes  of  Weather  and  Climate  in  Summer.  Hon.  Wea.  Rev. 

109.  In  press. 

Orianski,  I.,  1981:  The  Quasi-Hydrostatic  Approximation.  J.  Atmos.  Sci., 
l£,  572-582. 


19. 


Orlanski,  I.,  and  B.  B.  Ross,  1977:  The  circulation  associated  with  a  cold 
front.  Part  I:  Dry  Case.  J.  Atmos.  Sci.,  3^.,  1619-1633. 

Raymond,  W.  H.  and  H.  L.  Kuo,  1981:  A  Radiation  Boundary  Condition  for  Multi- 
Dimensional  Flows.  Submitting  to  J.  Comp.  Phys. 


Fig.  1 
Fig.  2 
Fig.  3 
Fig.  4 
Fig.  5 
Fig.  6 

Fig.  7 

Fig.  8 

Fig.  9 


Legend  of  Figures 

.  Wj  distribution  given  by  la  at  t  ■  15  hr. 

.  Same  as  Fig.  1  but  given  by  lb  at  t  ■  15  hr. 

.  Same  as  Fig.  1  but  given  by  2a  at  t  ■  15  hr. 

.  U j  distribution  given  by  2a  at  t  ■  51  hr. 

.  Wj  distribution  given  by  2b  at  t  -  15  hr. 

.  Vj  distribution  given  by  2b  at  t  ■  15  hr. 

colour  interval  4ms 
.  Wj  distribution  given  by  2b  at  t  *  24  hr. 
contour  interval  AW|  =*  10  cm  s 
dj  distribution  given  by  2b  at  t  ■  24  hr. 
contour  interval  AUj  -5ms 

.  Stream  function  i|>2  in  yz-plane  given  by  2b  at  t  -  24  hr. 

5  2  -1 

contour  interval  -  7.0  *  10  m  s  . 


I 


31. 


The  Dynamical  Structure  of  Squall-Line 
Type  Thunderstorms 


Abstract 

The  vertical  structure  of  the  squall-line  type  thunderstorm  is 
investigated  through  the  use  of  a  purely  two-dimensional  model  and 
the  flow  field  is  taken  as  occurring  in  an  unstably  stratified  atmosphere 
with  a  vertically  shearing  mean  wind.  A  case  study  of  a  real  squall-line 
appears  to  justify  this  simple  approach.  The  mechanism  which  leads  to 
the  most  characteristic  upshear  slope  of  these  storms  is  attributed 
to  liquid  and  solid  water  loading  and  attempts  are  made  in  this  work  to 
verify  this  hypothesis  by  integrating  the  dynamic  model  from  a  given 
initial  state  to  see  whether  it  will  evolve  to  such  an  arrangement, 
and  also  by  analysis  of  the  influence  of  a  reasonably  chosen  but 
fixed  liquid  water  distribution  on  the  vorticity  field  in  the 
vicinity  of  the  updraft-downdraft  interface.  The  preliminary  results 
indicate  that  the  proposed  mechanism  is  working. 
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1.  Introduction 

This  work  constitutes  a  part  of  our  effort  to  study  the  dynamics  of  the 
squall-line  type  thunderstorm  and  it  deals  only  with  the  usually  observed  up- 
shear  slope  of  this  type  of  storms.  The  main  thrust  of  this  research  is  to 
determine  the  cause  of  the  upshear  slope  of  the  updraft-downdraft  couplet  in 
the  squall-line  type  thunderstorms  under  the  influence  of  the  vertically 
shearing  mean  flow,  and  especially  to  find  out  whether  this  arrangement  can  be 
attributed  to  the  downward  drag  or  load  of  the  liquid  and  solid  water  and 
evaporative  cooling  from  the  falling  precipitation,  which,  according  to  the 
reasoning  of  the  senior  author,  indicates  that  the  upshear  arrangement  is 
most  efficient  energetically  for  the  disturbance  in  the  presence  of  precipi¬ 
tation  because  it  relieves  the  loading  from  the  updraft  and  turns  the  precipi¬ 
tation  drag  into  an  accelerating  force  for  the  downdraft,  whereas  with  either 
vertical  or  downwind  slope  these  forces  tend  to  destroy  the  updraft.  Further, 
it  is  evident  a  non-vertical  slope  can  not  be  maintained  without  a  vertically 
shearing  mean  wind  in  the  convective  storm. 

in  order  to  illustrate  the  nature  of  the  actual  squall-line  type  storms, 
a  brief  case  study  of  a  squall-line  which  moved  through  Chicago  area  recently 
is  presented  in  section  2.  This  is  followed  by  a  presentation  of  the  model 
and  the  numerical  method  used  to  solve  it  in  sections  3  and  4,  while  section  5 
discusses  the  results  that  have  been  obtained  from  the  model.  The  influence 
of  the  liquid  water  content  on  the  slope  is  discussed  further  in  section  6 

rv  n 

on  the  basis  of  the  vorticity  n  ■  /3x  -  .  Plans  for  further  improvement 

of  the  results  and  future  works  are  discussed  briefly  in  section  7.  while  a  few 
concluding  remarks  are  added  in  section  8. 
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2.  Squall  line  case  study 

In  the  late  evening  of  6  August  1 98 1  a  squall  line  developed  and  moved 
across  the  area  of  southern  Wisconsin  and  northern  Illinois.  Although  no 
severe  weather  was  reported  in  conjunction  with  this  squall  line,  it  did 
produce  high  winds  and  heavy  rain.  The  line  formed  along  a  trough  line  which 
ran  southwest  from  a  surface  low  in  northern  Wisconsin  to  northern  Mossouri. 

The  0000  GMT  7  August  1981  souding  from  Peoria,  IL,  is  shown  in  figure  1. 

Note  that  the  atmosphere  is  nearly  dry  adiabatic  below  900  mb  and  very  dry 
above  900  mb  except  for  a  few  thin  moist  layers.  This  sounding  shows  the 
atmosphere  to  be  "ripe"  for  convection,  and  with  a  wind  shear  of  about  40/m  sec 
through  the  depth  of  the  troposphere,  we  would  expect  the  convection  to  be 
organized. 

Figure  2  shows  the  Marsailles,  IL  radar  PPI  display  for  four  times  during 
the  squall  line  at  approximately  1  h  10  m  intervals,  starting  at  0037  GMT, 

7  August  1981.  At  the  earliest  time  shown,  there  is  a  line  of  broken  indivi¬ 
dual  cells  of  moderate  intensity.  As  the  cells  increased  in  intensity,  they 
became  elongated  and  began  merging  into  an  unbroken  line.  By  0148  GMT,  portion 
A  of  the  line  is  over  25  miles  long  but  less  than  10  miles  wide,  and  showed  a 
very  uniform  reflectivity  along  its  length.  The  storm  reached  its  peak 
intensity  at  approximately  0256  GMT.  At  this  time  the  line  was  very  long  and 
uniform.  The  line  remained  like  this  for  about  40  minutes  before  breaking  up 
into  the  weaker  ceils  shown  in  the  last  picture  at  0405  GMT.  In  all,  portion 
A  of  this  squall  line  was  visible  as  long,  thin,  coherent  line  of  at  least 
moderate  Intensity  for  over  2  hours. 

Figure  3  shows  the  hodograph  plotted  from  the  wind  data  of  the  0000  GMT 
Peoria  sounding.  The  speed  and  direction  of  portion  A,  and  the  orientation 


of  the  squall  line  where  determined  from  the  radar  pictures  and  are  also  shown 
on  this  diagram.  As  can  be  seen,  the  winds  were  somewhat  complicated  but 
show  a  generally  unidirectional  shear.  That  is,  there  Is  no  systematic 
veering  or  backing  of  the  winds  with  height.  One  also  sees  that  the  cells 
in  the  squall  line  move  in  a  direction  of  approximately  the  mean  wind. 

it  is  common  practice  in  studying  squall  lines  to  take  a  two-dimensional 
cross-section  perpendicular  to  the  squall  line.  This  "slab  symmetry"  approach 
assumes  that  the  line  is  uniform  along  its  length  so  any  motion  into  or  out 
of  this  plane  is  not  important.  In  this  case,  however,  the  line  is  oriented 
at  an  angle  of  nearly  45°  to  the  plane  of  the  mean  wind.  Thus,  a  cross-section 
taken  perpendicular  to  the  line  would  have  winds  into  and  out  of  the  plane  of 
the  same  magnitude  as  the  winds  in  the  plane.  Even  though  this  case  shows  a 
striking  uniformity  along  the  squall  line,  it  is  not  totally  uniform  and  so 
this  slab  symmetry  approach  does  not  seem  appropriate. 

A  schematic  representation  of  what  appears  to  be  happening  in  this  case 
is  shown  in  figure  4.  In  the  squall  line  the  updraft  and  downdraft  branches 
of  the  circulation  are  confined  to  the  plane  of  the  mean  wind,  which  is 
parallel  to  the  direction  of  motion.  The  line,  however,  is  oriented  at  an 
angle  to  this  plane  to  form  a  sort  of  "snowplow"  effect.  With  this  view  in 
mind,  it  Is  clear  that  the  plane  parallel  to  the  direction  of  motion  is  the 
physically  relevant  plane  of  cross-section.  We  shall  use  the  term  two-dimensional 
in  this  sense,  that  is,  the  circulation  is  approximately  confined  in  a  plane, 
and  therefore  refrain  from  using  the  term  slab  symmetry  since  it  does  not 
necessarily  imply  the  same  meaning. 

Taking  the  cross-section  parallel  to  the  direction  of  cell  motion  and 
defining  a  new  coordinate  system  with  x'  parallel  to  the  plane  and  y'  per¬ 
pendicular  to  it,  one  can  decompose  the  winds  into  the  components  u'  and  v' 
as  in  figure  5>  This  shows  the  v*  component  to  be  nearly  zero  In  the  mean 
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and,  therefore,  probably  Insignificant  to  the  organization  of  the  storm.  The 
u'  component  shows  the  storm  to  be  imbedded  in  an  environmental  shear  of  about 
3.3x10  ^  Sec  * . 

Ludlum  Cl 980)  defines  a  Richardson  number  for  storm  circulations  as 


-Ri 


CAPE 

(AU)2 


(2.1) 


where  CAPE  is  the  Convective  Available  Potential  Energy  and  can  be  expressed 
as  the  positive  area  of  a  lifted  parcel  on  a  skew-T  log-P  thermodynamic  dia¬ 
gram.  The  quantity  AU  is  the  algebraic  difference  of  the  horizontal  windspeed 
of  the  parcel  on  outflow  and  inflow,  or  more  simply,  the  difference  between 
the  winds  at  the  storm  top  and  the  surface.  Using  the  Peoria  sounding,  the 
Richardson  number  for  the  parcel  lifted  from  the  surface  is  -Ri  *  1 .98.  If 
the  surface  to  950  mb  average  values  of  temperature  and  moisture  are  used  we 
obtain  -Ri  *  1.45.  Ludlum  ( 1 980)  argues  that  -  Ri  <_  5  implies  organized  super¬ 
cell  type  convection  since  then  the  energy  due  to  shear  is  comparable  to  the 
energy  released  by  latent  heat. 

The  above  case  exhibits  the  features  which  are  indicative  of  the  squall¬ 
line  type  supercell  which  is  being  investigated  in  this  research  effort.  It 
is  quasi-steady,  maintaining  itself  with  little  change  for  a  period  of  2  hours. 
The  circulation  is  two-dimensional  in  the  sense  that  the  environmental  wind  is 
basically  confined  to  a  plane  parallel  to  the  direction  of  notion.  And,  the 
circulation  is  embedded  in  an  atmosphere  of  moderate  to  strong  shear,  yielding 
a  small  Richardson  number  for  the  storm. 

3.  Equations  and  boundary  conditions 

Instead  of  treating  it  as  an  initial  value  problem,  here  we  shall  seek  the 
time  dependent  model  and  allowing  an  initial  overturning  circulation  to  converge 
towards  a  steady  state.  This  is  an  established  method  in  this  type  of  problem 
(Roache,  1976).  There  is  also  no  boundary  between  the  updraft  and 
the  down-draft.  If  a  free-slip  boundary  was  present  between 
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the  updraft  and  downdraft,  it  would  represent  a  sheet  of  infinite  vorticity 
and  infinite  temperature  gradient.  Clearly  this  is  not  the  case  in  a  real 
thunderstorm  and  the  generation  of  vorticity  at  this  interface  by  the  finite 
temperature  gradient  is  a  major  dynamical  feature  of  the  storm.  As  well  as 
leading  to  a  more  realistic  model,  the  removal  of  the  interface  simplifies 
the  numerical  treatment  of  the  problem. 

The  basic  form  of  the  equations  is  the  same  as  that  given  in  Seitter 
(1980)  but  with  density  stratification  included.  They  are 
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where 

u,  w  *»  horizontal  and  vertical  velocities 
n  **  momentum  vorticity 
T  «  temperature 
q  “  water  vapor  mixing  ratio 
M  ■  liquid  water  content 

and 

p  *  Ps  exp(-z/HQ)  -  density 
T  ■  (1  ♦  0.6lq)T  ■  virtual  temperature 
F  *  dry  adiabatic  lapse  rate 
L  ■  latent  heat  of  condensation 
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Cp  ■  specific  Mat  at  constant  pressure  for  dry  air 

B(z)  -  exp  (z/H  .R/c  ) 

o  p 

R  ■  gas  constant  for  dry  air 
Hq  *  density  scale  height 

G  “  liquid  water  generation  term  (condensation  rate  for  q) 

V  ■  terminal  velocity  of  raindrops 

The  velocities,  u  and  w,  and  the  vorticity,  n>  can  be  written  in  terms 
of  a  streamfunctlon 


pu 
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n  - 
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where  V2  ■  -  +  -  .  The  liquid  water  generation  term  is  given  by 
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where  E  is  the  evaporation  rate,  a  constant,  and  qg  is  the  saturation  mixing 
ratio,  determined  from  the  Clausius-Clapeyron  equation.  The  terminal  fall- 
speed  of  raindrops  is  taken  as  a  constant  with  a  value  of  7  m/sec  (Seitter, 
1980). 


(3.6) 


A  schematic  of  the  model  circulation  is  shown  in  figure  6.  The  inflow 
characteristics  of  T,  q,  and  M  are  specified  for  both  the  updraft  inflow, 

0  <  z  <  zAu,  and  the  downdraft  Inflow,  <  z  <  H.  A  condition  on  n  is  also 
needed  and  is  taken  as  ■  0.  On  outflow  all  quantities  must  satisfy 


3A 
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A  ■  \p,  T,  q,  M,  and  n 


(3.7) 


The  value  of  ip  is  fixed  on  the  top  and  bottom  boundary  with  il> _ *  ip.  * 

top  bottom 

The  value  of  is  determined  such  that  i|/  ■  0  at  the  steering  level  z  *  z* 
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for  a  given  inflow  shear,  The  values  of  T  and  q  are  fixed  on  the  top  and  bottom 

except  when  precipitation  is  present.  Then  if  precipitation  has  reached  the 

lower  boundary,  a  method  similar  to  Schlesinger  (1973)  will  be  employed  to 

allow  the  boundary  to  cool  and  become  saturated  in  response  to  the  presence  of 

rainwater.  At  the  top  boundary  no  liquid  water  can  be  present  soH*0 

3M 

there,  on  the  bottom  boundary  we  set  =  0  to  allow  rainwater  to  "fall  through" 
the  surface. 

With  the  advective  terms  written  in  the  flux  form  no  boundary  value  for 
the  vorticity  need  be  specified  on  the  horizontal  boundaries  since  w  =  0  on 
them.  However,  if  the  Arakawa  Jacobian  is  used  or  if  a  diffusion  term  of  the 
form  vV2n  is  added  to  (3.1),  a  value  of  ri  on  the  boundary  is  needed.  Since 
this  model  is  intended  to  simulate  a  nearly  inviscid  atmosphere,  a  free-slip 
condition  is  appropriate.  Roache  (1976)  has  stated  that  the  proper  free-slip 
condition  on  the  vorticity  is  *  0  for  a  finite  difference  formulation. 

This  form  for  the  boundary  condition  gives  a  less  viscous  effect  of  the  boun¬ 
daries  than  the  normally  used  n  ■  0  condition,  even  though  the  latter  is 
analytically  correct. 

4.  Numerical  Method 

Equations  (3-1)  through  (3.4)  have  been  nond i mens i ona 1 i zed  and  written  in 
finite  difference  form,  with  the  velocities  u  and  w  written  in  terms  of  the 
streamfunction  using  (3- 5a) .  All  space  derivatives  were  written  in  centered 
difference  form.  Care  must  be  taken  to  write  products  of  quatities  in  a 
finite  difference  form  which  is  consistent  with  the  form  of  the  advective  terms. 
These  four  equations  and  ( 3 - 5b )  form  a  complete  set  in  q,  T,  q,  M  and  ip.  The 
method  of  solution  consists  of  starting  with  an  initial  field  for  each  of  the 
variables  and  integrating  (3>1)  through  (3.4)  forward  in  time  by  an  appropriate 
method.  Equation  (3.5b)  is  solved  by  successive-over-relaxation  (SOR)  for  the 
new  streamfunction  field  after  each  time  step.  This  process  is  continued  for 
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the  desired  time  or  until  a  steady  state  is  reached. 

Three  methods  of  time  integration  have  been  used  in  the  development  and 
testing  of  this  model.  A  forward-in-ttme  centered- in-space  method  was  used 
initially  as  a  simple  method  to  test  the  model.  This  method  was  soon  abandoned, 
however,  since  it  required  a  large  diffusion  term  for  linear  stability.  The 
inviscid  equations  were  also  integrated  using  the  two-step  Lax-Wendroff  method 
(Richtmeyer,  1962).  Although  this  method  has  been  used  with  success  in 
meteorological  problems  (Houghten,  et.  al.,  1966),  Lilly  (1965)  warns  that 
it  can  lead  to  severe  inaccuracies  in  certain  types  of  subsonic  fluid  flow 
problems.  This  seemed  to  be  the  case  with  the  current  problem,  as  it  was 
not  possible  to  eliminate  small,  but  accumulating,  errors  which  are  inherent 
in  the  Lax-Wendroff  method  and  which  eventually  destroyed  the  accuracy  of 
the  integration.  This  method  was  ultimately  abandoned,  although  some  impor¬ 
tant  results  were  obtained  while  it  was  being  used  (see  Section  6).  The 
very  accurate  method  of  Adams  and  Bashforth  (Lilly,  1965)  has  been  used  with 
success  both  in  integrating  the  truly  inviscid  equations  and  in  solving  the 
equations  with  an  added  viscosity  term. 

The  Adams-Bashforth  method  can  be  written  by  considering  the  equation 
for  some  quantity  A, 
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=  f 
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where  f  is  a  function  of  the  other  variables  and  their  derivatives.  We  let 
f^  be  the  value  of  f  evaluated  by  using  the  values  of  the  variables  at  time 
t  ■  nAt.  Then  the  value  of  A  at  time  t  ■  (n  +  1)  At  is  given  by 


A  ("■*•!) 


,(n) 


:(n) 


_  1  f(n-l) 
2 


] 


(4.2) 


This  method  is  slightly  amplifying  and,  therefore,  not  stable  in  the  strict 
sense  of  the  word.  The  amplification  factor  is  small,  however,  and  does  not 


40 


pose  a  problem  as  long  as  the  integration  is  not  extended  for  too  long. 

The  method  can  be  stabilized  by  adding  a  viscous  term  of  the  form 
vV2A.  This  is  included  by  writing  (4.2)  as  (Roache,  1976) 

>+l)  _  A(n)  A  f(n)  _  1  f(n-1)j  +  AtvV2^) 


1  f("-l 

2  T 


(4.3) 


Note  that  this  is  not  equivalent  to  including  a  viscous  term  in  f. 

Since  the  Adams-Bashforth  method  can  be  used  successfully  on  the  fully 
inviscid  equations,  the  viscosity  which  is  added  can  be  made  whatever  size 
is  appropriate  for  the  problem.  This  is  preferable  to  the  Lax-Wendroff  method 
that  has  a  built  in  numerical  diffusion  which  selectively  damps  small  wave¬ 
lengths. 

The  current  version  of  the  model  has  21  grid  points  in  the  horizontal 
and  11  in  the  vertical.  The  grid  spacing  is  even  and  taken  as  1  km  in  both 
the  horizontal  and  vertical,  giving  a  domain  of  20  km  by  10  km.  The  choice 
of  a  iO  km  depth  for  the  domain  was  somewhat  arbitrary.  Excluding  the  over¬ 
shooting  top  of  the  updraft,  observed  thunderstorm  circulations  extend  from 
the  surface  to  a  height  of  10  km  to  14  km  (Chisholm  1973,  Ludlum,  1980). 

The  depth  of  the  model  makes  very  little  difference  when  the  density  is  held 
constant,  but  in  simulations  in  which  the  variation  of  density  with  height  is 
included,  the  depth  of  the  model  may  have  to  be  adjusted. 

5.  Numerical  results 

As  discussed  in  Mon^.rieff  (1978),  the  case  of  a  neutrally  stratified, 
incompressible  atmosphere  with  M  =>  0  is  quite  simple.  In  this  case  the  vorticity 
equation  reduces  to 


d_ 

dt 


n  *  0 
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So,  vorticity  is  conserved  along  a  streamline.  If  we  consider  an  environment 
of  linear  shear,  then  on  inflow  the  vorticity  is  given  by 

n  ■  ■  constant  (5-2) 


since  ^  ■  0  on  inflow.  Thus,  the  vorticity  will  be  constant  throughout  the 
circulation,  and  the  steering  level  height  will  be  at  zAu  =  z^d  *  H/2. 

Moncrieff  shows  that  in  this  case  the  Interface  between  the  updraft  and  down- 
draft  will  be  vertical.  This  situation  Is  shown  in  figure  7* 

Figure  7  was  determined  by  performing  a  relaxation  on  a  vorticity  field 
which  had  constant  vorticity  everywhere  except  on  the  dividing  streamline. 

On  the  dividing  stream  line  the  vorticity  was  adjusted  to  yield  a  streamfunc- 
tion  value  equal  to  that  on  the  upper  and  lower  boundaries.  The  nondimensional 
streamfunction  values  marked  on  the  figure  correspond  to  a  linear  shear 
2  x  10  ^  sec  ' . 

The  inviscid,  incompressible  form  of  the  vorticity  equation  was  integrated 
forward  in  time  as  a  test  of  the  numerical  model  for  this  case.  The  flux  form 
of  finite  difference  was  used  for  the  advective  derivatives  and  the  Adams-Bash- 
forth  method  was  used  for  the  time  integration.  The  vorticity  field  of  figure  7 
was  used  as  the  initial  condition.  After  20  minutes  of  simulated  time  the  change 
in  this  field  was  less  than  the  initial  errors  In  it.  This  demonstrates  that 
the  flux  form  for  the  advective  derivatives  was  able  to  accurately  maintain 
this  known  solution  to  the  equations,  and  that  the  slight  amplification  inherent 
in  the  Adams-Bashforth  method  was  too  small  to  pose  a  problem. 

This  is  the  only  case  in  which  the  Adams-Bashforth  method  was  used  without 
diffusion  term  as  described  in  the  preceding  section.  Without  diffusion  of 
some  sort,  all  disturbances  must  be  advected  out  of  the  domain  and  very  small 
temperature  variations  can  lead  to  large  variations  in  the  vorticity. 
b.  Dry,  convectively  unstable  case 

In  order  to  further  test  the  model  by  reproducing  Moncrieff' s  results  of 


downshear  orientation,  the  equations  were  reduced  to  the  case  of  a  dry,  con 
vectively  unstable,  incompressible  atmosphere.  Then  the  equations  are 
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Since  r  is  taken  as  a  constant  equal  to  6.5  °C/Km,  a  constant  environmental 
lapse  rate  which  is  greater  than  this  yields  a  convectively  unstable  atmos¬ 
phere  and  exactly  models  the  equations  used  by  Moncrieff  (1978). 

Moncrieff  defines  a  type  of  Richardson  number  given  by 


(5.3) 


(5.4) 


(5-5) 


where  uq  is  the  horizontal  velocity  at  the  surface  in  the  updraft  inflow. 

This  is  related  to,  but  larger  in  absolute  value  than,  the  Richardson  number 
defined  by  (2.1).  Figures  8(a)  and  8(b)  show  the  results  of  Moncrieff1 s 
calculations  for  R  ■  H  using  his  analytical  and  numerical  models,  respectively. 
The  method  of  solution  for  the  analytical  result  has  been  described  previously 
(Seitter,  1 980) .  The  numerical  solution  was  found  using  a  method  similar  to 
the  current  model  except  that  a  staggered  grid  Lax-Wendroff  integration  scheme 
was  used. 

The  current  model  was  run  with  an  unstable  lapse  rate  which  gave  R  «  *s. 

The  initial  condition  was  given  by  the  circulation  shown  in  figure  7  and  an 

initially  adiabatic  atmosphere.  Diffusion  was  included  in  the  model  with  the 

2 

value  of  the  eddy  vis;osity  set  at  450  m  /sec.  While  this  is  about  an  order  of 
magnitude  larger  than  the  value  in  cumulus  clouds  (Tag,  1979) >  it  does  not  seem 
to  be  too  large  in  light  of  the  measurements  of  kinetic  energy  dissapation  in 
thunderstorms  (Frisch  and  Strauch,  1975),  and  is  much  smaller  than  the  numeri¬ 
cal  diffusion  included  in  many  previous  models  (Schlesinger,  1973).  It  was 
found  necessary  to  include  a  diffusion  term  because  without  it  the  generation 
of  negative  vorticity  at  the  updraft-downdraft  interface  by  the  temperature 
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gradient  there  formed  a  countered culat ion  at  the  interface. 

The  results  of  this  model  after  21  minutes  simulated  time  are  shown  in 
figure  9*  There  are  some  distinct  differences  between  this  result  and  those 
of  Moncrieff.  The  most  notable  difference  being  that  the  current  result 
shows  the  circulation  sloping  less  steeply  downshear.  Also,  the  steering 
level  heights  for  both  the  updraft  and  downdraft  are  located  at  z  =  H/2.  It 
is  felt  that  both  of  these  differences  are  a  result  of  the  inclusion  of  diffu¬ 
sion  in  the  model.  The  diffusion  of  temperature  weakened  the  gradients  and 
hence  the  production  of  vorticity,  and  the  diffusion  of  vorticity  "slowed  down" 
the  flow.  It  is  felt  that  this  circulation  represents  the  case  in  which  the 
instability  is  great  enough  to  overcome  the  viscosity  of  the  fluid  but  not 
great  enough  to  accelerate  the  flow  appreciably  and  thus  raise  the  steering 
level  height  zAu-  Thus,  the  modeled  flow  represents  a  flow  which  would 
correspond  more  directly  to  a  smaller  value  of  R  and,  hence,  would  slope  down- 
shear  less  steeply  (Moncrieff,  I 978) . 

Although  the  Lax-Wendroff  method  has  a  built  in  diffusion,  it  is  felt  that 
this  did  not  alter  Moncrieff's  results  in  the  way  it  did  in  the  current  model 
because  the  diffusion  in  the  Lax-Wendroff  method  is  highly  selective  toward 
short  wavelengths.  This  allows  the  method  to  prohibit  the  development  of  a 
counter-circulation  at  the  interface  while  allowing  the  rest  of  the  flow  to 
remain  virtually  inviscid.  This  kind  of  selective  damping  can  be  incorporated 
into  the  present  model  for  a  quantity  A  by  replacing  vV2A  with  K|V2A|V2A 
(Cullen,  1976).  It  is  felt  at  this  time,  however,  that  the  simpler  diffusion 
term  is  preferable. 

The  shape  of  the  interface  in  the  results  of  the  current  model  more  closely 
resembles  the  shape  of  Moncrieff's  analytical  result  (figure  8(a))  than  does 
Moncrieff's  numerical  result  (figure  8(b)).  This  is  felt  to  be  significant  and 
is  a  favorable  indication  as  to  the  accuracy  of  the  current  model. 


Although  no  details  are  given  in  either  his  paper  (Moncrieff,  1978)  or 
his  thesis  (Moncrieff,  1970)*  it  is  assumed  that  Moncrieff  must  have  used  a 
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much  smaller  grid  spacing  than  the  current  model  to  overcome  the  problems 
associated  with  the  Lax-Wendroff  method.  The  author  has  found  that  a  very  small 
grid  spacing  does  improve  the  Lax-Wendroff  method,  but  the  cost  of  increasing 
the  number  of  grid  points  substantially,  as  well  as  the  inherent  high  cost  of 
this  method  (over  twice  as  time  consuming  per  time  step  as  the  Adams-Bashforth 
method),  make  the  Adams-Bashforth  method  more  desirable. 

6.  The  mechanism  for  upshear  slope  in  the  light  of  the  vorticity  n 

it  has  long  been  recognized  that  the  upshear  slope  of  the  updraft  in  thun¬ 
derstorms  is  beneficial  to  their  maintenance  from  a  thermodynamic  point  of  view. 
However,  no  adequate  theory  has  been  proposed  which  explains  how  the  updraft  is 
able  to  oppose  the  shear  in  which  it  is  imbedded.  The  framework  of  the  present 
model  allows  a  dynamical  theory  for  the  upshear  slope  based  on  vorticity  argu¬ 
ments. 

Consider  the  incompressible  vorticity  equation 


3n  .  an  .  an  _  i  3Tv  .  3m 

?ttu3;  +  “5J  9Cr-  w  *  9  sr 
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(6.1) 


New,  consider  a  flow  pattern  as  shown  in  figure  10(a),  with  the  updraft-downdraft 
interface  vertical.  Let  the  updraft  have  a  positive  temperature  excess  and  the 
downdraft  have  a  negative  temperature  excess,  and  let  there  be  precipitation  in 
the  updraft,  in  figure  10(a),  the  precipitation  is  shown  in  the  lowest  half  of 
the  updraft  because  it  would  tend  to  be  heaviest  there  and  thus  have  the  most 
dynamical  effect  in  this  region.  As  can  be  seen  from  (6.1),  it  is  the  horizon¬ 
tal  derivatives  of  the  temperature  and  liquid  water  which  are  important.  The 
circulation  is  driven  by  the  temperature  distribution,  which  produces  negative 
vorticity  generation  (-|2.  <  o)  along  the  interface  and  hence  maintains  the  strong 
shear  across  the  interface.  Positive  vorticity  is  generated  (-g^-  >  0)  in  both 
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the  updraft  and  downdraft  branches  away  from  the  interface,  and  this  also  acts 

to  maintain  the  circulation.  Now,  consider  the  effect  of  the  liquid  water.  On 

dM 

the  downdraft  side  of  the  precipitation  column  (^  >  0),  positive  vorticity  will 
be  generated,  and  on  the  updraft  side  of  the  column  (^  <  0) ,  negative  vorticity 
will  be  generated.  There  are  two  ways  that  the  updraft  branch  can  respond  to 
this  negative  vorticity  generation:  1)  the  shear  can  be  reduced  to  make  the 
shear  vorticity  less  positive;  2)  the  streamlines  can  curve  less  sharply  as 
they  enter  the  updraft  to  reduce  the  curvature  vorticity.  Both  of  these 
responses  probably  take  place,  and  both  tend  to  cause  the  streamlines  to  enter 
the  updraft  at  a  more  upshear,  rather  than  vertical,  angle.  A  similar  argument 
leads  to  the  streamlines  in  the  lower  downdraft  becoming  more  sharply  curved  as 
they  exit  the  downdraft,  which  is  consistent  with  the  entire  circulation  sloping 
upshear. 

Consider  now  the  situation  shown  in  figure  10(b),  which  more  closely  resem¬ 
bles  the  flow  in  a  real  thunderstorm.  This  diagram  is  very  similar  to  the 
schematic  given  by  Ludlum  (1963)*  The  flow  is  still  driven  by  the  positive  and 
negative  temperature  excesses.  (Now,  the  downdraft  is  shown  to  have  a  negative  • 
temperature  excess  only  in  the  region  where  rain  water  can  cause  evaporation 
cooling,  as  is  the  case  in  real  storms.)  The  streamlines  show  less  curvature 
in  the  updraft  inflow  where  negative  vorticity  is  being  generated  and  more 
curvature  in  the  downdraft  outflow  where  positive  vorticity  is  being  generated. 
This  vorticity  generation  by  the  liquid  water  term  does  not  operate  unchecked, 
of  course,  it  opposes  the  effects  of  the  environmental  shear  which  tends  to  tilt 
the  circulation  downshear  until  a  balance  is  established. 

A  simple  means  of  testing  this  hypothesis  is  to  integrate  (6.1)  with  fixed 
distributions  of  temperature  and  liquid  water.  This  was  done  while  the  model 
was  in  the  Lax-Wendroff  phase  of  development. 

Figure  11  shows  the  results  of  the  integration  after  14  minutes  of  simulated 
time,  starting  with  an  initially  vertical  interface.  In  this  run,  a  temperature 


excess  of  —.25  #C  was  fixed  in  the  updraft  and  downdraft  regions  respectively 
and  no  liquid  water  was  present.  The  vertical  shear  on  inflow  was  2  x  10  ^ 
sec  This  should,  in  a  simple  way,  model  the  results  of  Moncrieff  (which 

were  modeled  more  accurately  in  the  previous  section).  Consistent  with  Moncrieff 
results,  the  flow  is  oriented  downshear.  The  flow  is  not  steady  at  this  time 
but  is  intensifying  and  beginning  to  form  a  closed  counter  circulation  between 
the  updraft  and  the  downdraft  (note  the  closure  of  the  96  streamline).  The 
steering  level  heights  Zjly  and  reflect  the  fact  that  the  flow  is  being 
accelerated  so  the  depth  of  the  inflow  must  be  greater  than  the  depth  of  the 
outflow. 

Using  the  same  temperature  distribution  and  the  same  initial  circulation, 
a  region  of  liquid  water  equalling  2g  m  ^  was  added  to  the  lower  half  of  the 
updraft.  Figure  12  shows  the  results  of  this  run  after  the  same  period  of 
time.  Now,  the  vorticity  generation  due  to  the  precipitation  has  had  a  marked 
effect  on  the  flow.  While  the  circulation  does  not  have  a  truly  upshear  orien¬ 
tation,  the  trend  is  definitely  to  oppose  the  effects  of  the  shear.  The  flow 
is  also  weakened  in  the  updraft  branch  by  the  effects  of  the  water  loading  and 
the  reverse  circulation  is  inhibited  except  in  the  region  above  the  precipita¬ 
tion. 

The  previous  results  seem  to  strongly  support  the  hypothesis  that  the 
upshear  orientation  of  the  updraft-downdraft  couplet  is  both  caused  and 
maintained  by  the  effects  of  the  liquid  water  which  is  produced  in  the  storm. 

This  would  explain  why  strong  thunderstorms  which  have  the  ability  to  produce 
large  amounts  of  rain  exhibit  an  upshear  slope  and  small  showers  and  cumulus 
clouds  slope  downshear. 

7.  Discussion  and  future  plans 

The  preceding  sections  have  shown  the  current  model  to  reproduce  the  results 
of  Moncrieff  fairly  well.  The  differences  can  be  explained  by  simple  reasoning 
based  on  the  differences  between  the  current  model  and  Moncrleff's  model.  The 
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simple  model  has  also  given  insight  into  the  mechanism  which  causes  the  upshear 
slope  of  the  updraft.  These  results  give  one  confidence  In  the  model  and  show 
that  the  proposed  extension  of  the  model  to  include  precipitation  is  warranted. 
This  section  will  discuss  the  improvements  which  have  been  made  in  the  model  and 
plans  for  its  use  in  further  investigations,  which  will  include  the  addition  of 
the  equations  for  water  vapor  and  liquid  water  to  the  system  of  equations  (3.1)- 
(3 - and  their  integration  to  a  steady  state.  To  improve  the  resolution,  the 
grid  spacing  will  be  cut  in  half.  It  is  felt  that  a  grid  spacing  of  no  greater 
than  500  m  in  both  the  vertical  and  horizontal  is  necessary  since  large  varia¬ 
tions  can  occur  over  the  distance  of  a  kilometer  or  two,  especially  in  the 
interface  region. 

It  is  hoped  that  the  compressibility  of  the  atmosphere  does  not  signifi¬ 
cantly  affect  the  qualitative  behavior  of  the  storm  structure.  This,  however, 
will  haveto  be  tested.  The  equations  were  proposed  in  the  anelastic  form  and  have 
been  finite  differenced  in  this  form  for  the  model.  Most  of  the  runs  are 
expected  to  be  carried  out  with  the  density  held  constant  due  to  the  relative 
ease  in  interpreting  the  results.  A  few  compressible  simulations  can  then  serve 
to  show  that  the  qualitative  structure  remains  the  same  as  well  as  providing  a 
more  realistic  picture  of  the  storm's  structure. 

The  moisture  equations  have  been  finite  differenced  and  tested  to  prepare 
for  their  inclusion  into  the  model.  Figure  13  shows  the  thermodynamic  diagram 
which  is  produced  by  the  equations  for  an  incompressible  atmosphere.  The  dry 
adiabatic  lapse  rate  is  taken  as  a  constant  with  the  value  9.76  °C/km.  The 
moist  adiabats  are  reproduced  quite  accurately  in  the  model.  Diagrams  like 
figure  13  will  allow  an  environmental  profile  to  be  plotted  and  the  corresponding 
Rl  to  be  calculated  for  each  case  using  (2.1)  and 
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where  6  is  the  potential  temperature  of  a  parcel  ascending  from  the  surface  and 

O 

0e  is  the  environmental  potential  temperature. 

The  situation  of  characterizing  a  flow  by  its  value  of  Ri  will  not  be  as 
easy  now  compared  to  Moncrieff's  model  in  which  constant  lapse  rates  for  both 
the  ascending  parcel  and  the  environment  were  assumed.  Now,  two  yery  different 
profiles  could  be  constructed  to  yield  the  same  numerical  value  of  Ri .  Differ¬ 
ing  profiles  of  environmental  moisture  can  also  be  expected  to  influence  the 
resulting  flow.  While  it  is  not  possible  in  this  more  complicated  model  to 
completely  characterize  a  flow  by  the  value  of  Ri,  this  still  is  expected  to  be 
an  important  parameter  and  will  be  continued  to  be  used  as  a  measure  of  the 
relative  importance  of  bouyancy  and  shear. 

The  first  case  to  be  run  with  the  complete  model  will  be  for  an  environment 
with  a  linear  temperature  profile  having  a  surface  temperature  of  300  K  and  a 
temperature  of  237  K  at  10  km  (y  ■  6.3  °C/km).  The  moisture  will  be  given  by 
90%  relative  humidity  throughout  the  atmosphere,  and  the  environmental  shear 
will  be  constant  with  a  value  of  4  x  10  J  sec  (AU  *  40  m/sec).  This  gives 
a  value  of  -Ri  =  1.0.  This  situation  is  not  meant  to  model  an  observed  thunder¬ 
storm  environment,  but  is  a  simple  environment  with  the  appropriate  Ri.  This 
case  should  provide  a  definitive  test  of  the  upshear  slope  hypothesis  and  provide 
a  benchmark  for  comparisons  of  runs  with  different  environmental  conditions. 

8.  Conclusions 

After  a  good  deal  of  effort,  a  numerical  method  which  is  reliable  and 
accurate  has  been  developed  and  the  basic  results  of  earlier  work  (Moncrieff, 
1978)  have  been  reproduced.  With  the  approval  of  a  grant  of  time  on  the  NCAR 
system,  it  is  now  possible  to  improve  and  complete  the  model.  The  full  model 
is  currently  being  readied,  and  after  testing,  it  will  be  run  on  its  first 
case.  This  case  should  confirm  the  hypothesis  which  was  presented  here;  that 
the  effects  of  the  liquid  water  act  to  produce  and  maintain  the  upshear  slope 
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of  the  updraft  which  is  thermodynamically  necessary  for  the  maintenance  of  the 

squall-line  type  thunderstorm. 
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Figure  2.  Marsailles,  IL.  Radar  PPI  display  for  four  times 
during  squall  line  passage.  The  reflectivity  contours 
are  for  levels  1,  3  and  5,  and  alternate  between  solid 
and  dotted  for  sequental  times.  Times  indicated  are  in 
GMT  on  7  August  1981.  The  distance  scale  along  the  90 
radius  is  in  mi les. 


Hodograph  of  the  envi ronmental  wind  taken  from  the  000 
ight  of  points  on  the  hodograph  is  shown  in  kilometers 
1  motion  and  orientation  of  the  line  is  also  shown. 


Figure  8  Results  of  Moncrieff's  analytical  (a)  and  numerical  (b 
model  for  R  =  -i.  Chain  line  in  (a)  represents  the  interface 
between  the  updraft  and  downdraft.  In  (b)  the  nondimensiona 
J*^e(’ature  e*cess  's  shown  as  a  chain  line  (From  Moncrieff, 


Figure  10.  Schematic  of  overturning  circulations.  The  solid 
lines  are  streamlines,  the  dashed  lines  are  isotherms  of 
positive  or  negative  temperature  excess,  and  the  cross- 
hatched  areas  represent  heavy  precipitation. 
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13-  Thermodynamic  diagram  produced  by  the  model  equations  for  an  incompressible 
tmosphere.  The  dashed  lines  are  dry  adiabats  and  the  solid  lines  are  moist 
diabats. 
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A  Radiation  Boundary  Condition  for  Multi-Dimensional  Flows 


Abstract 

A  'n'  dimensional  radiation  boundary  condition  for  application  at 
open  or  computational  boundaries  is  formulated  and  tested  on  three 
two-dimensional  problems  (n=2).  Two  of  these  problems  model  a  simple 
wave  propagation  and  possess  analytic  solutions  so  that  the  effective¬ 
ness  of  the  boundary  condition  can  be  measured  in  terms  of  a  RMS 
error.  A  more  subjective  analysis  must  be  used  in  the  final  problem 
which  is  the  the  simulation  of  an  atmospheric  cold  front.  The 
proposed  radiation  boundary  condition  requires  the  scalar  components  of 
the  phase  velocity.  A  formula  for  computing  these  components  is  given 
and  various  numerical  schemes  are  tested.  The  traditionally  used  one¬ 
dimensional  Sommerfeld  radiation  condition  is  recovered  when  n=l.  The 
higher  dimensional  radiation  boundary  condition  is  found  to  give  signifi¬ 
cant  improvement  over  the  one-dimensional  method  when  the  flow  is 


multidimensional . 


1.  INTRODUCTION 
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One  of  the  problems  facing  modelers  of  meso  and  other  small  scale  atmos¬ 
pheric  phenomena  is  that  in  these  finite  area  simulations  there  is  a  difficulty 
in  prescribing  lateral  boundary  conditions  since  no  true  physical  boundary  exists. 
The  nature  of  the  environment  outside  the  region  under  investigation  is  also 
unknown.  This  problem  has  been  forced  on  many  mesoscale  investigators  and  a 
variety  of  techniques  have  been  utilized  to  help  eradicate  this  difficulty,  e.g., 

(a)  a  Sommerfeld  radiation  condition; 

(b)  an  absorbing  boundary; 

(c)  one  sided  differencing  of  the  equations; 

(d)  various  other  types  of  extrapolation. 

Commonly,  these  procedures  are  utilized  on  some  very  complicated  problems  where 
analytical  solutions  do  not  exist  hence  the  impact  of  any  one  of  these  boundary 
condi tons  is  not  known  fully-  for  example,  Clark  [2]  using  different  expressions  for 
the  phase  velocity  associated  with  method  (a)  has  found,  for  flow  over  a  bell- 
shaped  mountain,  significant  variations  in  the  interior  calculations. 

What  is  needed  in  problems  where  advection  or  wave  motions  dominates,  as 
pointed  out  by  Orlanski  [9],  is  an  'open*  boundary  condition.  Such  a  condition 
entails  determining  if  the  'flow  pattern'  is  entering  or  exiting  across  a 
boundary.  In  the  latter  case  the  disturbance  should  be  allowed  to  propagate  out 
without  reflection.  It  Is  in  this  spirit  that  Orlanski  [9]  and  Pearson  [11] 
proposed  to  use  the  following  form  of  the  Sommerfeld  radiation  condition  at  the 
boundary: 


(l) 


where  $  is  the  variable,  C  the  phase  velocity,  t  the  time  and  nr  the  coordinate  per¬ 
pendicular  to  the  boundary  in  question.  Pearson  [11]  proposed  estimating  C  from 
a  linearized  dispersion  relation  while  Orlanski  [9]  proposed  to  determine  C 
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locally  and  hence  to  predict  the  boundary  value  of  <{>  without  finding  the  disper¬ 
sion  relation  which  as  a  rule  is  unknown.  It  is  clear  that  with  the  Orlanski 
approach  Eq.  (1)  can  predict  accurately  one  dimensional  motion  but  it  does  not 
seem  adequate  to  represent  the  higher  dimensional  flows.  In  this  study  a  *n' 
dimensional  radiation  condition  is  proposed  and  tested  together  vlth  different 
techniques  in  evaluating  the  components  of  the  phase  velocity  on  three  two-dimen¬ 
sional  problems  two  of  which  possess  analytic  solutions.  For  more  information 
on  other  boundary  techniques  see,  e.g.,  Engquist  and  Majda  [3],  Gustafsson  and 
Kreiss  [4],  Rudy  and  Strikwerda  [13]  and  Schubert  et  al.  [14], 

Typically  in  limited  area  or  mesoscale  atmosphere  studies  the  boundaries  are 
placed  as  far  as  possible  from  the  center  of  activity.  To  test  and  evaluate  the 
proposed  boundary  condition  a  somewhat  different  philosophy  is  taken  in  that  we 
study  the  distortion  as  a  phenomena  nears  and  passes  through  a  boundary. 


2.  A  'n'  DIMENSIONAL  RADIATION  BOUNDARY  CONDITION 
According  to  the  definition  of  the  phase  velocity  £  of  the  field  4>,  in  general 
we  have 


||  -  -  -  -  E  c,  Uj*' .  i  -  l,2,...,n. 


(2) 


where  C.  is  the  component  of  ?  in  the  direction  of  Xj.  For  our  three-dimensional 
space  we  have  1  <  n  <  3.  We  shall  use  this  relation  as  our  general  radiation 
condition  at  the  boundary  provided  t  is  directed  outward  from  the  boundary,  and 
determine  C.  by  applying  the  governing  equation  <{>  which  we  shall  write  formally  as 

||  "  -  F(xj ^ - ,xn>t,4»).  (3) 

Thus,  on  equating  the  two  expressions  of  given  by  (2)  and  (3)  we  obtain 

E  c*  H-  -  FU,  ••  •»*  .  <*> 

jTj  dxj  •>  n» 

For  the  one-dimensional  problem,  Cj  can  be  obtained  directly  from  this  equation 


provided  F  is  known.  However,  for  the  n  dimensional  problem  with  n  >  1,  we  need 


n  equations  to  determine  all  C.  uniquely.  To  reduce  the  number  of  unknowns  to  one, 
we  make  use  of  the  property  of  C  implied  by  the  relation  (2),  namely  £  is  parallel 
to  V<J>-  In  accordance  with  this  property  we  assume  that  we  have  the  optimal  relation 


a 
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where  the  constant  of  proportionality  factor  a  is  the  same  for  all  i.  On  substitu¬ 
ting  this  relation  in  (4)  we  then  obtain  at  «=  F/.^  and  therefore 


C!  " 


(6) 


Equivalently,  we  may  write  that  each  component  of  CxV$>=(?  is  satisfied  and  together 
with  Eq.  (4)  this  yields  n  equations  and  n  unknowns  with  (6)  as  the  solution.  Thus, 
the  value  of  every  Cj  at  the  boundary  can  be  determined  from  the  adequate  forms 
of  and  F,  such  as  equating  F  to  the  value  of  "34>/gt  at  one  grid  point 

inside  the  boundary  and  one  time  step  earlier,  as  proposed  by  Orlanski  [9].  This 
approach  will  be  tested  together  with  the  techniques  in  determining  F,  by  applying 
it  first  to  two  two-dimensional  problems  where  solutions  are  known  analytically, 
and  then  to  the  more  complicated  two-dimensional  cold  front  problem  without  known 
solution. 


3.  THE  THREE  TESTING  PROBLEMS 


? rob  1 em  A 

A  two-dimensional  advection  problem  In  non-dimensional  form  given  by 

0<x<1  ,  0<y<! ,  t  >  0  ,  (7) 

is  solved  for  u  *»  v  ■  1  with  initial  conditions 

u(x,y,o)  *  AsIn27rxsin2Try.  (8) 

The  boundary  conditions  along  the  inflow  boundaries  at  y*0  and  x*0  are  made  to 
satisfy 

u(x,y,t)  -  Asin2tr(x-t)s in2” (y-t)  ,  (9) 

which  also  represents  an  analytic  solution  of  (7)  for  every  point  of  the  domain. 


69. 


This  solution  projects  a  pattern  which  moves  at  a  45°  angle  toward  the  upper 
right  corner. 

Eq.  (7)  is  solved  numerically  using  a  leap-frog  scheme,  with  Robert's 
filter,  for  a  time  step  of  .01  on  a  21  by  21  equally  spaced  grid  having 
Ax=Ay=.05,  and  applying  the  open  boundary  conditions  of  the  form  discussed  in 
section  2  at  the  outflow  boundaries  along  x=l  and  y=l.  The  constant  A 
is  assigned  the  value  of  100.  The  RMS  error  is  computed  in  the  traditional 
manner. 

Problem  B_ 

Rossby  waves,  a  commonly  known  large  scale  atmospheric  motion  feature 
represented  by  the  solution  of  the  barotropic  nondivergent  vorticity  equation 
whose  linearized  version  is  given  by 


3r'  -  3r'  -3r*  ,  3f 

“St  +  U  JT  +  V  3y  "  0  » 

where  u  and  v  ere  the  basic  currents  in  x  and  y  directions,  f  is  the  Coriolis 
parameter  and  u' ,  v' ,  and  are  the  perturbation  velocities  and  perturbation 
vorti  city  which  are  expressed  in  terms  of  the  perturbation  stream  function  ij;  by 

„•  V'  3  v  - 1^-  -  ~  o<x<x 

u  3y,  v  ^  ,  5  7T  3y  ^  °-X-XN  ’ 


(10) 


°<y<yM  .  (11) 

It  is  well  understood  that  Rossby  waves  owe  their  existence  to  the  variation  of 

the  Coriolis  parameter  f  with  latitude. 

For  our  testing  purposes,  i.e.,  calculations  limited  to  a  finite  region, 

we  will  assume  a  solution  of  the  form 

i|>  **  A  cos  k  (x  -  C  t)  cos  m  (y-C  t)  ,  (12) 

a  y 

C*  *  "  A(k2+m2)  cos  k(x-C  t)  cos  m(y-C  t) ,  (13) 

x  y 

where  k  •  2tt,.  is  the  zonal  wave  number,  m  “  2ir.,  is  the  meridional  wave  number 
'  x  '  y 

while  L  and  L  are  the  wave  lengths  in  the  x  and  y  directions,  respectively.  The 
x  y 

x  component  of  the  phase  velocity  is 


here  0  =  3f/3y 


Cx  =  u  -  (k2+m2) , 

,  while  the  y  component 


is  g i ven  by 
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The  initial  conditions  are  given  by  Eqs.  (12)  and  (13)  for  t  =  0  while  the 
boundaries  along  x  =  0  and  y  “  0  are  also  obtained  from  the  analytical  solution. 

The  known  outflow  boundaries  along  x^  =  3000  km  and  y^  =  2250  km,  for  both  the 
vorticity  and  streamfunction  calculations,  are  used  for  testing.  The  unknown 
streamfunction  boundary  values  are  computed  using  the  radiation  formulation  and 
the  phase  velocity  components  obtained  from  the  vorticity  equation.  A  leap  frog 
scheme  is  used  to  solve  Eq.  (10)  and  S0R  technique  is  used  for  Eq.  (11).  A  21  by  16 
evenly  spaced  grid  is  used  in  the  x  and  y  directions  respectively  for  a  step 

size  of  150  km  along  with  a  time  step  of  2000  s.  In  addition  0  =  l.6Xl0-,Ism 

g  2  -  **1”  “1 

A  =  10  m  s  ,  u  =  15  m  s  ,  v  =  0  or  5  m  s  ,  L  »  6000  km  and  L  *  3000  km. 

x  y 

Problem  £ 

For  the  final  problem  we  choose  a  purely  two-dimensional  anelastic  moist  cold 
front  model  to  simulate  the  circulation  associated  with  an  atmospheric  cold  front 
(Ross  and  Orlanski  [12]).  A  cold  front  represents  a  propagating  disturbance  that 
can  not  be  described  completely  as  wave  motion  thus  problem  C  differs  from  problems 
A  and  B.  The  governing  equations  for  the  (x,z)  plane  are  of  the  form 


3n 

at 


Jf'J'.CM'l)  + 


s_ 

e 

o 


36 

v 

w 


-  f- 


3v 

3z 


1_(  jRy 

3z  '  <)z ' 


3c 

3x 


06) 


If  -  -  J(<M0n)  ♦  ozvw  +  fiu-ug)  +  |j  (vk  If)  .  ), 


If  -  -  J(*.ooe)  *  ♦  f  +  |;(KeK||)  +  fj(Ke|j)  +  . 


e  v  3u 


36v  .  L6 


It  “  ‘  J(*’aoq)  +  °zqw  +  l7(KeK^}  +  ll(Ke  I?  “  6  » 


-  J(i|>,aoc)  +  (K@K  ||)  +  (K#  ||)  +  6  +  o  we  . 


07) 

08) 

09) 


3c 

7t 


(20) 


71 


3w 

Here  u,  v  and  w  are  the  velocity  components  in  x,  y  and  z-dl rections,  n  “  w/3x 
-  *s  the  y-component  of  vorticity.ip  is  the  momentum  stream  function  in  xz- 

plane,  8  is  the  potential  temperature,  g  is  gravity  acceleration,  q  and  c  are  the 
water  vapor  and  cloud  water  mixing  ratios,  0v  ■  0(1  +  . 608q),  is  the  basic  state 
geostrophic  wind,  L  is  the  latent  heat  of  condensation,  c,,  the  specific  heat  at 
constant  pressure.  6  is  the  condensation  rate,  v  and  Ke  are  the  eddy  viscosity 
and  eddy  diffusivity  coefficients,  K  is  a  constant  used  to  enhance  the  horizontal 
diffusivity,  and  cr  a  -  3£npo/3z  is  the  stratification  factor  of  the  undisturbed 
density  p^.  Details  of  how  to  calculate  the  condensation  rate  6  can  be  found  In 
Ross  and  Orlanski  (12)  and  Kuo  and  Qian  (7). 

In  terms  of  the  stream  function  ip,  the  velocities  u  and  w  and  the  vor- 
ticity  n  are  given  by 


3(J>  dip 

p  T,.2iS  ♦»!*♦<,  |t  , 

o  -  2  «2  Z  02 


(21a, b) 
(22) 


3z 


where  a  “1,  and  p  =  p  ,  exp(-o  z)  is  the  vertically  varying  density, 
o  /PQ  o  surf  2 

In  our  calculations  the  grid  spacing  is  such  that  Ax»Az.  Under  these  condi¬ 
tions  the  streamfuncti'on  has  traditionally  been  calculated  by  neglecting  the 
32^^x2  term  in  Eq.  (22)  yielding  the  well  known  hydrostatic  approximation.  Orlan¬ 
ski  [10]  has  introduced  a  quasi  hydrostatic  approximation  in  which  a  correction 
containing  part  of  the  non-hydrostatic  contribution  is  added  to  the  hydrostatic 
solution.  In  the  Orlanski  [10]  procedure  the  streamfunction  is  represented  as 


m 


j  «o  ■* 

and  the  components  of  i}>  are  obtained  from  the  following  equations: 


(23) 


32i|> 

3z! 


dip 

2  *  °z  sr  -  - 


(2*1.) 
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at  z  =>  H,  iJ>o  «■  ¥  , 


-  0,  j  >  1 


r*i 
2 


+  0  fii  -  -  .  j  >  * .  (24b) 

32-  2  ^  3x2 

Here  i|>o  is  the  hydrostatic  component  and  the  Vj's,  j  >  1,  are  an  approximation 
to  the  non-hydrostatic  contribution.  The  associated  boundary  conditions  are: 

at  z  *  0,  ifij  -  0,  for  all  j  ,  (25a) 

(25b) 
(25c) 

The  series  for  converges  rapidly  provided  Ax  >  Az.  Taking  m  -  2  Is  sufficient 
for  our  purposes.  Thus  the  solution  of  (22)  is  obtained  successively  from  Eqs. 

(24a, b). 

At  the  lateral  boundaries  the  right  hand  side  of  Eq.  (24b)  must  be  approxi- 

3iJ'5 

mated.  Instead,  as  an  alternative,  we  choose  to  let  “  0  ,  j  :>  1 ,  l.e.,  the 
non-hydrostatic  terms  do  not  contribute  to  the  vertical  velocity  at  the  lateral 
boundary.  Otherwise,  there  are  no  lateral  boundary  conditions  required  by  this 
procedure. 

For  the  problem  in  general  at  the  lower  boundary,  z  ■  0,  slip  boundary 

30 

conditions  are  utilized,  e.g.,  T}  =  0,  =  0  and  v  satisfies  the  thermal  wind 

relation  given  by 


3v  g  30 

SI  "  le0JI  * 


(26) 


At  the  upper  boundary,  i.e.,  H  »  14  km.,  we  have  a  rigid  lid  so  that  w  »  0,  thus 
i  is  a  constant,  and  in  addition  n>8  and  v  keep  their  initial  vertical  gradients. 
Conditions  at  the  lateral  boundaries  are  left  for  experimentation. 

Initially,  the  thermal  wind  relation  (26)  is  taken  as  satisfied  by  v  and  0, 
andUgis  given  some  vertical  variation.  Generally,  our  initial  conditions  and  the 
values  used  in  the  eddy  viscosity  formulation  as  well  as  all  pertinent  details  for 
q  and  c  are  similar  to  those  used  by  Ross  and  Orlanski  [12]  and  thus  will  not  be 
presented  here. 


The  general  numerical  approach  is  composed  of  a  lumped  finite  element 
scheme  with  a  leap  frog  time  integration.  The  Arakawa  representation  of  the 
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Jacobian  J(ij>,&),  where 

f-Hif  <! 

is  obtained  when  bi  linear  el  enents  (chapeau)  are  used  in  the  finite  element 
formulation.  This  was  first  identified  by  Jesperson  [5].  The  diffusion  terms 
are  lagged  one  time  step  and  a  Robert's  filter  is  used  to  reduce  the  tendency 
of  time  splitting.  A  grid  utilizing  Ax  **  20  km  and  Az  ■  500  m  is  used  along  with 
a  100  s  time  step.  The  value  of  ¥  used  in  Eq.  (25b)  depends  only  on  the  initial 
u  velocity  field,  the  latter  is  a  function  of  the  z  coordinate  only. 


4.  NUMERICAL  PROCEDURE 

Radiation  boundary  conditions  of  the  form  of  Eq.  (2),  for  the  case  n  ■  2  are 
to  be  used  in  the  model  problems  at  the  appropriate  boundaries  as  Indicated  in 
section  3.  For  notational  purposes  rewrite  Eq.  (2)  for  the  n  »  2  case  as 

+  c  +  c  |i  -  o  .  (2 

3t  x  3x  y  3y 

Eq.  (28)  can  be  used  to  predict  the  value  of  <J>  on  an  open  boundary  provided  the 

value  of  C  and  C  are  known.  These  values  can  be  obtained  via  Eq.  (6),  e.g., 

X  V  -1 


The  relationship  between  C  and  C  is  given  by 

x  y 

Cx  "  Cy  3*/3x  /a*/3y  .  i 

Two  methods  of  evaluating  Eqs.  (29)  and  (30)  are  apparent  and  are  now  presented. 
An  Extrapolation  Approach 

The  values  of  4>  at  any  lateral  boundary,  say  j=$(xN,yj)  ,j*2, . . .  ,M-1 ,  as 
predicted  by  Eq.  (28)  must  be  determined  numerically.  If  Eq.  (  28)  is  evaluated 
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by  a  leap  frog  approach  for  a  time  step  centered  at  t  ■  T  then  the  correct  formu¬ 


lation  for  C  Tand  C  Tshould  correspond  to  the  phase  velocity  components  at  time 

jk  y 

t  and  centered  correctly  in  space.  Provided  these  C  's  and  C  's  are  not  rapidly 

x  y 

varying  in  time  and  space  an  approximation  based  on  previous  interior  values  of 
<{>  can  be  made,  e.g.,  Orlanskt  [9]  used  cjj  ■  c|J_J  in  the  one  dimensional  radiation 
boundary  condition  Cn  *  1),  where  ■  C(^,t»T).  Other  schemes  are  possible  but 
this  approach  has  the  advantage  that  it  is  independent  of  the  numerical  procedure 
used  in  the  interior,  i .e. ,  to  solve  Eq.  C 3) •  However,  various  improvements  can 
be  made  and  these  will  be  indicated  later. 

Assuming,  like  Orlanski  [9],  that  F  Is  given  by  centered  at  T  -1 

and  K  "1  we  find  that  at  the  right  most  boundary  (>^  ,y^  ) ,  excluding  corners,  of  the 
rectangle  region  bounded  by  o<x<xN  and  o£y£yH  that  Eqs.  (29)  and  (31)  can  be 
approximated  by 


cx  *  *  rst^-ij '  C2i,j]li 


c  -  c  -P 
y  x  ^y 


M  ■ 


where 


[^4 


(32a) 


(32b) 


(32c) 


3y  "  [*N-i.j+i  ■  *N-!.j-i]/2Ay 


(32d) 


Thus  each  of  the  C  's  and  C  's  are  located  spatially,  in  the  upstream  sense,  at 
x  y 

(XN_1  yj^  an<*  a^out  T  “  1  time  step.  The  phase  velocity  components  can  be 
T+l 

used  to  predict  .  via,  e.g.,  an  implicit  formulation  of  Eq.  (28),  yielding. 


■  T+l  .  0  -  '.“'W  , 


'N-I.j 


-  ht.j*i ;  ’itj-i  _ 


where  D  ■  I  +  CxAt/Ax  and  with  the  restriction  given  in  Orlanski  (9)  expanded  so 


that 


Cx“JCx  ’  °-Cx-  ^/At  . 
'  ^7  At  ’  Cx-  ^/At  , 


(34a) 


Ay/At  *  Cy  >  Ay/At  * 


V  K  ,-^<Cy<Ay/At  ^ 


-  Ay/At  *  Cy  <  "  Ay/At  * 

Similar  formula  can  be  written  down  for  the  lateral  boundaries  at  (x.  ,yu), 

■  n 

etc.  If  corners  are  to  be  computed  then  both  C  and  C  must  be  computed  using 

x  y 

one-sided  differencing.  Note  that  the  Orlanski  [9)  formulation  for  the  one 
dimensional  radiation  condition  is  recovered  above  if  Cy  **  0.  Also,  if  the 
properties  at  an  inflow  boundary  are  known  then  Eq.  (33)  need  not  be  used  in  that 
si tuation. 

Bannon  [1]  and  Miller  and  Thorpe  [8)  have  both  proposed  an  upstream  time 
differencing  scheme  to  evaluate  the  one-dimensional  Sommerfeld  radiation  boundary 
condition.  The  latter  authors  suggest  using 


where  r  *  C  At/A  ,  o  <  r  <  1 ,  and  r  is  determined  from 


(<-\  ■  -  Q  . 


In  addition,  Miller  and  Thorpe  [8]  performed  a  truncation  error  analysis  and  found 
that  Improved  accuracy  is  obtained  when 


r  ■  (ci  ■  -  c.)  *(f„  -  ■  c) 

-  (Ci  -  C!)/(*5-2  -  Cl)  • 

We  will  make  comparisons  between  both  the  leap  frog  and  upstream  time  differencing 


schemes. 


Note  that  Eqs.  (36)  and  (37)  require  an  interior  value  evaluated  for  time  step  t+1 


/b 


If  the  Interior  calculations  are  made  using  a  leap  frog  or  some  explicit  scheme 
then  the  computed  solutions  at  time  t+1  can  be  utilized  in  the  lateral  boundary 
calculations.  Under  these  conditions  t  may  be  replaced  with  T+1  in  Eqs  (32a) 
through  (33).  The  phase  velocity  components  are  then  centered  correctly  in 
tire  for  the  leap  frog  integration  scheme  given  above.  For  some  problems  it 
is  possible  to  compute  the  phase  velocity  correctly  in  both  time  and  space. 

Equating  To  The  Equation  Technique 

Again  assume  a  leap  frog  time  integration  scheme  for  Eq.  (28).  If  in  Eqs 
(2S)  and  (30)  It  is  possible  to  evaluate  F(x,y,t,$),  described  by  Eq  (3),  at 
the  boundary  at  time  T  using  one-sided  differencing  then  and  are  obtained 
centered  spatially  at  (x^_^  yj)  and  at  the  t  time  step.  If  F  is  a  complica¬ 
ted  function  then  this  one-sided  differencing  approach  will  be  more  involved  than 
the  extrapolation  procedure.  In  addition  one  sided  differencing  may  itself 
introduce  large  errors  especially  if  the  equation  contains  terms  that  are  in  a 
state  of  near  quasi  balance,  e.g.  the  geostrophic  balance  condition  in  the 
pri-nitive  equations  of  motion  and  the  thermal  wind  balance  condition  in  the 
vorticity  equations.  This  error  can  be  removed  to  a  certain  extent  provided 
hicner  order  finite  differencing  approximations  are  used  at  the  boundaries.  It  is 
best  to  avoid  this  type  of  error  if  possible* Nevertheless  in  our  model  problem  A 
•  and  B  this  technique  can  be  used  to  help  measure  the  error  introduced  by  using 

the  extrapolation  procedure  described  above.  Thus  given  the  values  of  C  and  C 

T+  1 

for  Eqs  (29)  and  (30)  the  value  of  is  obtained  as  before  via  Eq  (33). 

5.  DISCUSSION  OF  RESULTS 

In  order  to  test  the  'n'  dimensional  radiation  boundary  condi t ion,  three  dis¬ 
tinct  two-dimensional  problems  have  been  proposed.  For  problems  A  and  B  the  use 
of  boundary  conditions  at  the  known  outflow  boundaries  over  specifies  these  problems 
since  they  are  first  order  in  each  of  the  independent  variables.  However  this  over 
specification  is  one  possible  method  that  can  be  used  to  test  the  ramifications  of 
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any  technique  that  might  be  proposed  for  use  at  an  open  or  computational  boundary. 


Res u Its  From  Problem  A 

Starting  with  problem  A  we  show  in  Table  1  the  RMS  error  after  one  hundred 
time  steps.  The  results  using  both  a  one  and  two  dime  *!onal  radiation  condition 
are  shown  for  two  methods  of  determining  the  phase  velocities.  The  extrapolation 
approach,  described  in  Eqs  (3ia-d) ,  as  indicated  goes  back  one  time  step  from  the 
last  known  value  at  the  boundary  and  Interior  to  the  boundary  one  space  step  to  deter¬ 
mine  the  C  's  and  C  's  while  the  equating  to  the  equation  technique  evaluates 
x  y 

these  phase  velocity  components  using  one  sided  finite  differencing  of  terms  in 
the  governing  equation.  Orlanski  19]  first  proposed  the  extrapolation  procedure 
with  the  one  dimensional  radiation  boundary  condition.  Under  the  extrapolation 
approach  note  that  there  is  a  five  fold  reduction  in  the  RMS  error  by  increasing 
the  radiation  boundary  condition  from  one  dimension  (n=l)  to  two  dimensions  (n*2) . 

The  RMS  error  is  further  reduced  for  both  the  one  and  two  dimensional  radiation 
conditions  by  using  the  equating  to  the  equation  technique.  For  this  problem 
the  two  dimensional  radiation  condition  has  exactly  the  same  form  as  the  equation 
being  solved  so  it  is  not  surprising  that  the  RMS  error  is  essentially  identical 
to  that  given  by  using  one  sided  finite  differencing.  Limiting  the  radiation 
condition  to  one  dimension  increases  the  error  substantially.  Using  the  known 
analytical  solution  on  all  boundaries  does  not  reduce  the  error,  as  shown  in 
Table  1,  due  to  inaccuracies  in  the  interior  numerical  calculations. 

In  Fig.  1  values  of  for  the  first  fifty  time  steps  are  shown  at  location 
,yJ2)*  The  solid  line  is  for  the  one  dimensional  radiation  condition  while  the 
dashed  curve  is  for  the  two  dimensional  case.  Both  are  for  the  extrapolation  leap 
frog  method  of  determining  the  C  's  or  C's.  The  dotted  curve  represents  the  C  's 

X  X 

obtained  using  the  two-dimensional  equated  to  the  equation  technique.  For  the  two- 
dimensional  cases  all  values  are  to  be  compared  against  the  analytically  predicted 
value  of  I. 
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Substituting  the  analytical  solution  into  the  one  dimensional  boundary  condi¬ 
tion,  Eq  (1),  yields  an  analytically  predicted  time  dependent  solution  for  the 
phase  velocity,  i.e.. 


I  + 


tan2Tr(x-t) 

tan2Tr(y-t) 


which  becomes  infinitely  large  in  magnitude  when 


(38) 


2iv(y-t)  •* +]it,  j  =  0,1,... 
or 

2ir(x-t)  »  +fvn  ,  h  »  1,3,5...  • 

2 

This  explains  why  in  the  numerical  procedure  the  value  of  C  (solid  line)  fluctuates 

rapidly  and  is  quasi  periodic  in  time.  Nevertheless  the  value  of  C  must  be  re- 

Ax 

stricted  (Orlanski ,  19] )  i.e.,  o  <  C  <_  otherwise  substitution  of  C  back  into 
the  formula  (Eq  33,  n~l)  to  predict  the  new  boundary  value  of  the  dependent  varia¬ 
ble  u  would  numerically  make  no  sense. 


Results  for  Problem  Bi 

In  Table  II  some  RMS  errors  are  given  for  problem  B,  the  Rossby  waves,  after 
a  total  of  one  hundred  time  steps  or  55-55  hrs.  The  results  are  presented  in  a 
format  similar  to  that  given  in  Tabie  I  except  two  cases  are  given.  The  first  two 
columns  of  the  RMS  errors  for  the  vorticity  and  streamfunction,  respectively,  are 
for  when  the  true  phase  velocity  is  entirely  in  the  x  direction,  i.e.  *  o, 
while  columns  three  and  four  are  for  when  the  phase  velocity  has  components  in 
both  coordinates.  The  latter  occurs  when  v  is  non  zero. 

In  Table  II  note  that  when  v  -  o  the  one  and  two  dimensional  methods  give 
approximately  the  same  RMS  errors  in  each  category  used  to  determine  the  phase 
velocity.  The  extrapolation  approach  gives  the  largest  errors,  as  compared  to 
the  equating  to  the  equation  technique,  and  predicts  changes  along  the  y*^ 
boundary  when  in  fact  no  changes  occur.  When  the  mean  flow  contains  components 
in  both  the  x  and  y  directions  the  two  dimensional  radiation  condition  again 
becomes  superior  as  seen  in  columns  three  and  four.  Also  notice  that  the  two- 
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dimensional  equating  to  the  equation  technique  and  one  sided  finite  differencing 

give  nearly  identical  solutions.  This  impltes  at  the  known  out  flow  boundaries 

that  the  components  of  the  phase  velocity  are  predicting  outflow  and  thus  satisfy 

the  restrictions  imposed  by  Eqs.  (34a, b). 

In  Fig.  2  values  of  along  the  boundary  x»x^  are  displayed  at  the  fiftieth 

time  step.  The  solid  and  dashed  curves  are  C's  for  n=I  and  C  's  for  the  n=>2  cases, 

x 

respectively,  obtained  by  the  extrapolation  approach.  The  dotted  curve  is  found 
from  the  two-dimensional  equated  to  the  equation  technique.  The  analytically  pre¬ 
dicted  value  for  the  components  of  the  phase  veloctty  are  C  =C  and  C  =C  . 

xx  y  y 

Substituting  the  analytical  solution  for  the  vorticity,  Eq  (13),  into  the  one 
dimensional  radiation  condition  gives  an  analytically  predicted  value  for  the  phase 
velocity,  i.e.. 


c  “  Cx  +  Cy  ^  tan  m  (y-Cyt^an  k(x-Cxt)  . 


(39) 


Here  our  previous  finding  is  again  repeated  since  this  equation  predicts  that  in 

the  one  dimensional  case  the  phase  velocity  can  take  on  values  much  in  excess  of 

•  Ax 

the  acceptable  upper  limit,  i.e.,  C  ■  and  values  much  less  than  the  lower 
limit  C  =  0.  The  latter  Is  obtained  inspite  of  the  fact  that  the  flow  is  con¬ 
tinuously  outward.  The  sharp  spike  in  the  solid  curve.  Fig  2,  reflects  this  large 
variability  and  is  commonly  found  at  least  at  one  grid  point  at  almost  every  time 
step  and  often  exceeds  the  maximum  magnitude  allowable.  It  is  clearly  seen  in 
Eq- (39)  and  in  Tables  I  and  II  that  theuseof  the  one  dimensional  radiation  boundary 
condition  incurs  greater  errors  as  the  flow  becomes  increasely  multi -dimensional . 

In  Table  III  the  RMS  errors  associated  with  the  upstream  time  differencing 
schemes  of  Miller  and  Thorpe  [8]  are  given.  The  simplier  of  their  two  schemes, 
computed  using  Eq  (36)  for  n*l ,  gives  nearly  identical  results  with  the  one-dimen¬ 
sional  (n»l)  case  evaluated  using  Eqs  (32a-d)  with  x  replaced  by  x  +  1 .  Utilizing 
the  interior  solutions  evaluated  at  time  x  +  l  In  the  lateral  boundary  calculations 
does  not  change  the  RMS  error  appreciatively  either  positively  or  negatively  as 


seen  by  comparing  Tables  II  and  HI.  The  interpolation  of  r  and  hence  C  by  Eq  (37)  jj 

ii 

does  show  definite  improvement  except  for  the  streamfunction  calculations  when 

'  i 

there  is  two-dimensional  mean  flow.  With  further  integrations  the  large  error  in  ] 

the  latter  category  will  eventually  deteriorate  the  vorticity  solutions.  j 

Some  severe  storm  modelers  have  had  success  using  a  constant  phase  velocity 

in  the  one-dimensional  radiation  condition,  e.g. ,  see  Klemp  and  Wilhelmson  [6]. 

Table  IV  shows  the  RMS  errors  when  C  and  C  are  each  held  fixed.  Using  the 

x  y 

analytically  predicted  values  of  C  ■  12.082  and  C  ■  0  or  5  in  Eq  (33)  gives 

X  7 

RMS  errors  essentially  identical  to  the  one-sided  finite  differencing.  From 
Table  IV  we  also  see  that  for  one  dimensional  flow  the  error  occurred  by  over 

l 

i 

estimating  C  is  less  than  when  C  is  underestimated.  However  this  pattern  is  j 

X  X  I 

I. 

not  clearly  reproduced  when  the  flow  is  two-dimensional.  With  the  higher  dimen¬ 
sional  radiation  condition  it  may  also  be  acceptable,  under  some  circumstances, 
to  use  a  fixed  or  constant  value  for  each  component  of  the  phase  velocity  provided 
enough  information  is  available  to  determine  the  nearly  correct  magnitudes  and 
di rect ions. 


Results  For  Problem  £ 

Problem  C  requires  the  numerical  solution  of  several  equations.  Even  though 
complicated  by  the  release  of  the  latent  heat  the  procedure  is  nevertheless  straight 
foreward  except  for  the  conditions  to  be  used  at  the  open  boundaries.  At  the  lateral  boun 
daries  we  tested  a  variety  of  different  conditions  for  the  various  equations.  These 
tests  show  that  the  procedure  used  by  Clark  [2],  and  many  others,  is  very  satisfac¬ 
tory,  i.e.,  the  velocity  component  normal  to  the  lateral  boundary  Is  computed  at 
the  boundary  from  the  radiation  condition  while  all  other  dependent  variables  have 
zero  normal  gradients.  Thus  in  our  problem  at  the  lateral  boundaries  6,q,c  and  v 
satisfy  ^(  )  -  0  while  the  radiation  boundary  condition  is  applied  to  the  vorticity 
equation  from  which  the  normal  velocity  component  is  computed  via  the  streamfunc¬ 
tion. 

Our  results  for  the  atmospheric  cold  front  simulation  are  very  similar  to 
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those  obtained  by  Ross  and  OrlanskI  I 1 2} .  Details  of  the  various  fields,  e.g., 
vorticity,  streamf unction,  etc.,  are  very  involved  and  it  Is  not  easy  to  gauge  from 
these  the  influence  that  the  lateral  boundaries  might  play.  A  clearer  picture 
showing  the  Influence  of  the  various  radiation  boundary  condition  formulations 
is  gained  by  examining  plots  of  an  averaged  quantity  or  spacial  norm  verses  time. 
Figs.  3  and  4  display  our  choices. 

In  Fig.  3  we  plot  u'  verses  time  where  u*  Is  the  norm  of  the  perturbation 

velocity  as  defined  by 

N  M 

■(TiT  S  j5  |ui.j'U9I  •  m 

Here  u  is  the  computed  x  component  of  the  velocity  field  while  Ug  is  the  geostrophic 
wind  which  Is  a  function  of  the  vertical  coordinate  only.  This  norm  is  choosen 
because  it  allows  for  an  easy  interpretation  of  the  magnitude  of  the  perturbation 
velocity. 

In  Fig.  3  the  solid  curve  displays  the  results  obtained  using  the  two  dimen¬ 
sional  radiation  condition  (n=2)  for  the  extrapolation  techniques  while  the  dashed" 
dot  curve  is  the  Orlanski  procedure  (n*«l)  and  the  dashed  curve  Is  the  Mi  1 ler-Thorpe 
technique  using  the  improved  estimates  of  r  as  given  in  Eq  (37).  The  dotted  curve 
which  is  essentially  identical  to  the  solid  curve  uses  phase  velocity  components 
(n=2)  centered  correctly  at  time  t  for  the  leap  frog  scheme  but  still  computed  by 
the  extrapolation  procedure.  In  these  calculations  only  the  radiation  boundary 
condition  has  been  changed.  Also  the  area  of  computation  coincides  with  the 
region  over  which  the  averaging  Is  performed.  The  grid  Is  as  defined  in  section  3 
for  the  x  and  z  coordinate  system  utilizing  N®55  and  M»29.  The  above  mentioned 
curves  are  to  be  compared  and  contrasted  with  the  plus  (+)  curve  which  gives  an 
average  over  the  same  area  but  when  the  calculations  are  performed  on  a  larger 
area  having  N»73. 

Two  features  are  pronounced  in  comparing  the  various  curves.  First,  there  is 
some  adjustment  to  the  initial  conditions  and  all  the  curves  show  an  increase  in 
magnitude  which  reaches  a  maximum  at  sixteen  hours.  During  this  process  the  curves 
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are  very  similar.  After  the  adjustment  process  the  plus  curve  shows  that  a  quasi 
steady  state  exists.  Here  the  small  variations  are  due  to  the  presence  of  internal 
gravity  waves.  The  second  feature  isnow  apparent  since  the  other  curves  reveal 
that  the  lateral  boundary  has  a  large  impact  if  the  disturbance  is  too  close.  The 
Mi  1 ler-Thorpe  procedure  is  particularly  susceptible  since  boundary  influences 
cause  the  value  of  u1  to  more  than  double  the  maximum  magnitude  observed  in  the 
plus  curve.  The  two  dimensional  radiation  boundary  condition  (solid  or  dotted 
curves)  displays  the  smallest  increases  for  times  greater  than  thirty  hours.  The 
small  differences  between  the  curves  computed  using  the  one  and  two-dimensional 
radiation  boundary  conditions,  before  thirty  three  hours,  can  be  explained  by  the 
fact  that  away  from  the  frontal  zone  the  u  component  of  the  velocity  is  at  least 
two  orders  of  magnitude  larger  than  the  vertical  velocity  component.  Thus  the 
flow  is  essentially  one  dimensional  except  near  the  front.  This  is  clearly  seen 
by  the  magnitudes  of  w  in  Fig.  h. 


In  Fig.  k  values  of  w,  where 


N  M 

5-  J_  £  £ 

N*M  i=l  j=l 


w, 


U 


(41) 


are  shown  for  calculations  using  the  same  radiation  conditions  as  in  Fig.  3  except 
now  the  dashed  curve  is  the  Mi  1 ler-Thorpe  upstream  technique  computed  using  Eq(36). 
Again  the  lateral  boundary  effects  are  clearly  evident  after  thirty  three  hours  in 
the  leap  frog  schemes  which  give  very  similar  results  and  after  twenty  six  hours 
in  the  upstream  procedure  (dashed  curve).  The  latter  technique  generates  vertical 
velocity  averages  twice  as  large  as  the  leap  frog  scheme.  The  exact  reason  the 
upstream  procedure  is  less  effective  is  unclear  in  light  of  its  performance  in 
problem  B  (Table  III).  The  flow  pattern  is  however  much  more  complicated  in 
problem  C  and  it  seems  reasonable  that  a  boundary  procedure  that  uses  the  same 
time  integration  scheme  (leap  frog)  as  used  in  the  interior  calculations  would  be 
the  most  compatible  for  long  time  integrations. 
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In  Eqs.  (17)  through  (20)  w  appears  explicitly  in  one  term  in  each  equation. 

In  our  calculations  we  compute  w  from  calculated  streamfunction  values  via  the 

formula  w  *  aQ  *n  finite  element  representation. 

Along  the  lateral  boundaries  however  the  finite  element  technique  predicts 

each  w  using  streamfunction  values  in  a  one  sided  differencing  scheme  from  six 

grid  points.  To  test  whether  the  reduced  accuracy  associated  with  the  one  sided 

differencing  scheme  might  be  enhancing  the  error  at  the  lateral  boundaries  we 

gw 

replaced  the  direct  calculation  of  w  with  the  condition  that  -g — =  o.  This 

greatly  reduces  the  noise  generated  at  the  lateral  boundaries  as  seen  by  comparing 

3w 

values  of  w  in  Table  V,  computed  using  =  o,  against  our  earlier  calculations  of 
w  in  Fig.  k  in  which  w  was  calculated  at  the  lateral  boundaries.  All  categories 
of  radiation  boundary  conditions  show  improvement  with  now  almost  no  differences 
between  the  n=l  and  n=2  cases.  However  the  leap  frog  approach  still  remains 
superior  to  the  up-stream  technique. 

In  our  lumped  finite  element  scheme  every  term  in  Eqs.  (17)  through  (20)  is 
expressed  in  a  nine  point  configuration  except  for  the  time  term.  Hence  the 
contribution  from  every  term  that  includes  w  can  be  interpreted  in  terms  of  the 
streamfunction  as  representing  even  more  grid  points  since  the  calculation  of 
each  w  involves  \p  on  nine  grid  points.  Consequently  using  w  instead  of  the 
streamfunction  representation  explicitly  obviously  results  in  some  smoothing. 

This  fact  is  reflected  in  the  values  of  w  shown  in  Table  VI.  In  the  last  two 
lines  in  Table  VI  it  is  clear  that  when  the  calculations  are  performed  using 
ao  explicitly  in  place  of  w  the  values  of  w  are  nearly  a  magnitude  larger 
after  forty-eight  hours.  The  differences  are  small  during  the  first  six  hours 
but  as  latent  heat  is  released  at  grid  point  locations  in  the  condensation  process 
the  differences  grow.  Unfortunately  the  intensified  flow  characteristics  obtained 
when  using  the  streamfunction  representation  explicitly  also  enhances  ail  noise  and 


the  solutions  become  very  noisy  after  twenty  simulation  hours.  Thus  It  {5 
difficult  to  test  radiation  conditions  when  the  meteorological  fields  deter lo- 
rate  beyond  the  point  of  interpretation.  Nevertheless  the  result  for  w 
displayed  in  Table  VI  agrees  essentially  with  our  previous  findings. 

To  determine  if  the  error  at  the  lateral  boundaries  can  be  reduced  even 
further  we  tested  several  other  ideas.  For  example,  in  an  attempt  to  center 
the  phase  velocity  component  correctly  in  space  we  computed  two  columns  of 
components  adjacent  to  the  boundary  and  then  utilized  a  Taylor  series  expansion 
to  predict  the  value  at  the  boundary.  These  calculations  did  not  significantly 
improve  the  results  presented  above.  Using  higher  order  approximations  of  the 
derivatives  also  did  not  significantly  change  our  results.  Also  averaging  all 
quantities  over  six  grid  points,  as  computed  via  the  finite  element  method  withbi' 
linear  basis  elements,  yields  essentially  the  same  results.  As  previously 
mentioned  the  equating  to  the  equation  technique  cannot  be  utilized  In  this 
problem  because  of  the  internal  balance  between  two  derivative  terms  In  the 
vorticity  equation,  i.e. ,  the  thermal  wind  relation. 

6.  SUMMARY 

A  'n'  dimensional  radiation  boundary  condition  has  been  tested  on  three 
two  dimensional  problems.  When  the  flow  is  outward  across  a  lateral  boundary 
and  when  the  pattern  of  movement  is  multi -dimensional  the  proposed  radiation 
condition  has  been  found  to  be  clearly  superior  to  various  formulations  of  the 
traditionally  used  one  dimensional  Sommerfeld  radiation  condition.  For  outflow 
the  'n'  dimensional  radiation  boundary  condition,  as  proposed,  is  equivalent  to 
one-sided  differencing  of  the  governing  equation  provided  the  components  of  the 
phase  velocity  are  correctly  centered  in  time  and  space.  When  this  centering 
procedure  is  not  feasible  or  gives  unrealistic  results  extrapolation  procedures 
provide  an  alternative  technique  to  determine  the  phase  velocity  components. 
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TABLE  I 


RMS  errors  after  100  time  steps  for  problem  A. 

Boundary  Condition _ RMS  error 

(I)  Extrapolated 

|i+c||-=0  (la)  2.024 

T 

+  C  +  C  |i»0(lb)  .3956 

3t  x  3x  y  3y 

(II)  Eauated  to  the  equation 

4£  +  c|$  *0  (Ha)  .4224 

at  oHj. 

d?  +  r  +  C  a  n  (I lb)  too-s 

3t+cx33c  y3y  u  -1383 

(III)  One  sided  finite  .1380 

di fferences 

(IV)  Analytical  exact  value  .2020 
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TABLE  II 

RMS  errors  after  100  time  steps  for  problem  B. 


u  *  15  m  s'1 

u  =  15  m  s'1 

<i 

n 

o 

v  =  5  m  s'*1 

Boundary 

Condition 

Vortfcfty*  10  V1 

Streamfunction.  10® 

2  -1 
m  s 

Vorticity,  10'5 
s_1 

Stream- 
function, 
106  m2s 

(la) 

.7336 

.8650 

.5971 

1.5072 

(lb) 

.1308 

.7180 

.1521 

.2878 

(Ha) 

.0112 

.0296 

.0674 

.6995 

(Hb) 

.0112 

.0270 

.0475 

.1929 

(III) 

.0113 

.0270 

.0448 

.1977 

TABLE  III 

RMS  errors  after  100  time  steps  for  upstream  and  leap-frog 
radiation  schemes  used  in  problem  B. 


u  »  15 

ms’1 

u  *  15  m  s"1 

o 

n 

i> 

v  =  5  m  s'1 

Schemes 

Vorticity*  10“5s"1 

Streamfunctlon,  I06 

2  -1 
m  s 

Vorticity.  10”5 
s'1 

Stream 
function  , 
10®  m2  s”1 

Up  stream 

Eq.  (36) 

.128 

.8166 

.561 

1.4631 

Eq.  (37) 

.0298 

.1245 

.1300 

1.6006 

Leap  Frog 

n  =  1 

.1285 

.8171 

.5864 

1.5106 

n  *  2 

.1287 

.6861 

.1642 

.3212 
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TABLE  IV 

RMS  errors  after  100  time  steps  for  constant  phase  velocity  components 
In  the  radiation  condition  used  in  problem  B 


Boundary 

Condition 

u  =  15 

v  =  0 

m  s'1 

u  =  15  m  s'1 
v  =  5  ms'1 

cy 

Cx 

Vorticity* 

10-5s-1 

Streamfunction,  10® 
n£s~^ 

Vorticity,  10 '®s-1 

Stream 
f  uncti  on, 

10s  m2  s~T 

0.0 

12.082 

.0114 

.0298 

.6846 

1.2197 

5.0 

12.082 

.0434 

0.1974 

0.0 

5.0 

.1521 

1 

.6614 

.8748 

1.4602 

0.0  25.0 


.0640 


3658 


.6116 


1 .6686 
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TABLE  VI 

«  -2 

Values  of  w  (10  )  computed  for  problem  C  when  w  is  replaced  with  a  3^ 

0  oX 

in  Eqs  (17)  through  (20)  or  when  w  is  used  explicitly. 


1 

4 

Time  (hrc 

8 

») 

16 

24 

28 

48 

using  ip 

Up  stream 

.104 

.203 

.339 

2.431 

2.098 

2.565 

3.577 

£ 

• 

c r 

LU 

Eq.  (37) 

.104 

.203 

.341 

2.584 

2.353 

2.418 

1.489 

Leap  frog 

n  =  1 

.104 

.206 

.342 

2.374 

2.459 

1.897 

. 

1.502 

n  =  2 

.104 

.207 

.342 

2.376 

2.464 

1.929 

1.556 

Large  area 

CM 

ii 

c 

.104 

.202 

.305 

2.344 

1.918 

1.818 

1.491 

using  w 

Large  area 

n  *  2 

, 

.098 

.184 

i 

.222 

.452 

1 

.273  j 

.283 

.211 

List  of  Figures 
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Fig.  I 


Fig.  2 


Fig.  3 


.  Values  of  the  phase  velocity  components  are  displayed  for  the  first 
fifty  time  steps  for  location  ^ 1 2^  *n  Pr°biem  A.  Values  computed 
using  the  one  dimensional,  n=l ,  (solid  curve)  and  two  dimensional, 
n=2,  extrapolated  (dashed)  and  equated  to  the  equation  technique 
(dotted)  are  to  be  compared  against  0^=1. 

.  Same  as  Fig.  1  except  for  problem  B  and  computed  for  the  fiftienth 

time  step  along  the  right  lateral  boundary.  The  known  analytical 

solution  is  shown  by  the  C  •  C  curve. 

xx 

.  For  problem  C  values  of  u1  verses  time  are  shown  when  in  the  radiation 

condition  n**2  (solid  and  dotted  curves),  for  n=l  (dashed  dot)  and  for 

an  upstream  time  differencing  scheme  (dashed  curve)  using  Eq.  (37). 


Fig.  4.  Same  as  Ffg.  3  except  for  w  and  the  upstream  scheme  uses  Eq.  (3b). 
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Simulation  of  Laboratory  Vortex  Flow 
by  Axisymmetric  Similarity  Solutions 


Abstract 

A  similarity  approach  is  utilized  to  investigate  a  simple 
axisymmetric  steady-state  model  of  the  convergence  region  of 
a  laboratory  vortex.  The  resulting  simplified  set  of  equations 
are  solved  for  a  range  of  swirl  angles  by  varying  the  tangential 
or  radial  velocity  component  at  the  outer  rim.  By  increasing  the 
swirl  angle  the  flow  is  found  to  go  from  a  one  cell  to  a  two  cell 
configuration,  i.  e.,  the  vertical  velocity  changes  from  everywhere 
positive  to  negative  in  the  vicinity  of  the  axis.  Correspondingly 
the  vertical  vorticity  maximum  moves  from  the  axis  toward  the  radius 
of  maximum  tangential  velocity,  making  the  flow  barotropical ly 
'instable  with  respect  to  unsymmetric  perturbations. 
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1.  introduction 


A  variety  of  rotating  or  swirling  motions  are  observed  in  the  atmosphere. 

Of  particular  interest  dynamically  are  the  tornado  type  vortices.  Various 
theories  on  the  nature  of  the  internal  dynamics  within  these  strong  vortex  flows 
have  been  discussed  (e.g.,  see  Davies-Jones  and  Kessler,  197*0,  but  actual  mea¬ 
surements  are  almost  impossible  to  obtain  due  to  the  severity  of  the  motions.  To 
answer  some  of  the  major  questions  a  number  of  laboratory  vortex  experiments  have 
been  conducted,  e.g.,  Turner  and  Lilly  (1963),  Ying  and  Chang  (1970),  Ward  1972, 
Jischke  and  Parang  (197*0,  Leslie  (1977),  Church  et  al .  (1977),  etc.  One  of  the 
advantages  of  the  modern  vortex  generator  is  its  ability  to  reproduce  the  multiple 
vortex  phenomena  in  an  environment  where  it  can  be  intensely  studied,  e.g..  Ward 
(1970,  1972).  Fujita  (1971,  1972)  first  proposed  that  the  cyloidal  markings 
observed  in  tornado  damage  paths  are  produced  by  these  secondary  rotations.  The 
presence  of  multiple  vortices  within  some  tornadoes  has  now  been  widely  accepted 
and  analysed  in  various  studies,  e.g.  Agee  et  al .  (1975)  and  Forbes  (1978). 

Because  of  the  success  of  the  laboratory  simulations,  a  few  elaborate  numeri¬ 
cal  reproductions  of  the  vortex  generator  have  been  attempted,  e.g.,  Harlow  and 
Stein  (197**)  and  Rotunno  (1977,  1979).  These  studies  certainly  have  helped  to 
identify  various  features  of  the  dynamics,  but  questions  still  remain,  e.g.,  what 
is  the  relationship  between  the  vertical  motion  near  the  axis  and  the  inflow  angle 
at  a  large  distance  from  the  axis,  and  what  is  the  nature  of  the  instability  that 
allows  multiple  vortices  to  generate.  To  help  refine  answers  to  these  and  other 
questions,  a  simple  analytical-numerical  steady-state  model  is  investigated. 


l 
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2.  Formulation  of  the  Problem 

Our  goal  Is  to  simulate  In  a  simple  manner  a  steady-state  vortex 
similar  to  those  produced  In  modern  laboratory  vortex  generators  and 
to  study  its  behavior.  These  laboratory  vortices  are  created  mechanically 
in  a  specially  designed  arrangement  which  contains  an  exhaust  fan  and  a 
rotating  screen.  Air  is  drawn  radially  inward, by  the  fan  through  the 
rotating  screen  till  it  reaches  a  central  cylinder  where  it  is  allowed  to 
rise  through  the  cylinder  and  Is  expelled  at  the  top. 

The  rotating  screen  imparts  a  background  angular  momentum  which  is 
nearly  conserved  as  the  air  moves  radially  Inward  giving  rise  to  a  con¬ 
centrated  vortex  core.  The  laboratory  experiment  is  constructed  so  that  it 
is  geometrically  and  dynamically  similar  to  conditions  found  in  the  atmosphere 
during,  for  example,  a  tornadic  event.  This  similarity  is  established  by 
requiring  that  the  non-dimensional  numbers  that  govern  the  flow  in  the 
mechanically  generated  vortex  be  in  the  same  range  as  those  observed  in  the 
atmosphere.  Lewellen  (1962)  found  three  nondimensional  parameters  that  govern 
the  nature  of  the  swirling  flow.  These  are:  the  radial  Reynolds  number 
Rer  «  Q/2irv  , 
an  internal  aspect  ratio 
a  -  H/R  , 
and  a  swirl  ratio 

S  ■  tan8/2a  . 

Here  Q  is  the  volume  flow  rate  per  unit  axial  length,  v  is  the  viscosity 
coefficient,  H  the  inflow  depth,  R  the  radius  of  convergence  and  8  is  the 
swirl  or  inflow  angle.  It  is  known  that  the  swirl  ratio  and  radial 
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Reynolds  number  describe  the  dynamics  of  the  flow  while  the  aspect  ratio 
describes  the  geometry. 

If  all  three  of  these  numbers  are  of  equal  importance  then  the  geometry 
of  the  simulator  would  have  to  be  considered  and/or  duplicated  in  any 
numerical  simulation.  Davies-Jones  (1973)  however  points  out  that  the 
nature  of  the  swirling  flow  is  dominated  by  the  swirl  ratio.  Because  of 
this  there  is  some  hope  that  a  restricted  investigation  ignoring  the  geometry 
and  only  dealing  with  just  part  of  the  vortex  simulator,  e.g.  the  region 
where  convergence  takes  place,  might  be  successful  in  explaining  some  of 
the  general  behavior  that  is  observed  in  the  laboratory  vortex  and  thus 
consequently  in  the  atmosphere.  It  is  in  this  context  that  we  adopt  a 
similarity  approach  to  simplify  the  governing  equations.  The  procedure  has 
been  orchestrated  so  that  we  can  examine  the  vortex  flow  configuration, 
as  described  in  the  simplified  similarity  equations,  for  a  range  of  swirl 
angles . 


3.  The  Mathematical  Development 

Assuming  that  the  motion  of  the  homogeneous  fluid  is  steady  and  axl- 
symmetrlc  and  that  the  Coriolis  force  can  be  neglected  in  comparison  with 
the  centrifugal  force,  then  the  swirling  flew  satisfies  the  following 
equations  of  motion  in  cylindrical  coordinates  (Kuo,  1966,  1969): 


“<|-r 


a  v 


0) 


jd..iu  +  v7|u 

3  r  3z  r  p  3  r  v  I 


3  w 


1  3 


(2) 


(3) 
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where 


3z 

r  Is  the  radial  coordinate,  z  the  vertical  coordinate  positive  in  the  direction 
opposite  to  the  gravitational  acceleration  represented  by  g,  and  u,  v  and  w  are 

the  radial,  the  tangential  and  the  axial  component  of  velocity,  p  is  the  pressure 
and  p  is  the  density. 

The  velocity  components  u  and  w  can  be  expressed  in  terms  of  the  stream 
function  defined  by 

u  •  "  ?fi  •  w  "  f  •  (5a, b) 


From  Eqs.  (2)  and  (3)  we  obtain  the  following  equation  for  the  azimuthal  vorticity 


(6a) 


(6b) 


Equations  (6b)  and  (1)  together  form  a  closed  system  for  the  dependent  variables 
v  and  <|)  (after  using  (5a, b)  and  (6a)  to  eliminate  u,  w  and  c) . 
a.  A  Similarity  Approach 

Following  the  spirit  of  the  similarity  approaches  utilized  by  Kuo  (1966, 
1969)  the  variables  are  first  non-dimensional  I  zed  by  the  following  transformations: 
1(1  ■  U#r  h$*  ,  z  -  hz* 

,  v  ■  Usr#v*  and  x  ■  r2Mr2 


m  ■  U  r  m* 

9  9 


(7) 
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Here  r$  is  the  radius  of  the  rotating  screen  at  the  outer  edge  or  rim,  h  is 
a  measure  of  the  effective  depth,  rQ  is  some  measure  in  the  radial  direction 
and  U$  is  the  radial  velocity  at  the  rim  while  m  ■  vr  denotes  the  angular 
momentum.  Note  that  v*,  ip*,  m*,  z*  and  x  are  all  dimensionless. 

The  dimensionless  equivalent  of  Eq.  (l)  when  written  in  terms  of  m* 


becomes 


3<P*  3m*  .  3ip*  3m*  «  .  r  d2m*  X-1  32m*| 

-  jFr-3irt3ir^-2''*i'^r+6  • 


while  Eq.  (6)  reduces  to 


m*  |$  -  6x2  [iv*  ^  -  I?  *  Ig.  (S ht) 


m  3z*  r  x  3x  '^c^z*  3z*  v  x*x|  ’  w/ 

where  the  subscripts  z*  and  x  denote  partial  differentiations  and 

6  -  h2r  -2  ,  02**  .  x  £ti  ♦  J-' 

°  3x2  3z*2 

and  D4^*  »  D2(D2ip*).  Observe  from  (9)  that  the  meridional  flow  is  influenced  by 
m*  only  through  the  variation  of  m*  with  z*.  To  include  this  effect,  we  follow 
again  the  basic  premise  established  by  Kuo  (1966,1969)  by  expanding  the  flow 
variables  in  power  series  expansions  in  z*and  including  the  zeroth  order  terms 
mQ  and  uq  in  m*  and  u*.  It  can  readily  be  shown  that  only  the  even  order  terms 
of  z*  in  m*  and  u*  will  contribute  to  these  two  variables  and  therefore  we  write 
m*  and  ip*  In  the  following  forms 

m*  m  v*[m0  +  6z*2mj  +  (Sz*2)2!^  +...],  (1C 

**  .  2v*z*(Fq  +  6z*2F,  +  (6z*2)2F2  +  ...J  .  (11 

Both  m*  and  tp*  will  be  well  defined  provided  the  higher  order  terms  are  decreasing 
in  importance. 

From  Eqs.  (5a, b),  (7),  ( 10)  and  (II)  we  arrive  at  the  following 


t 


— — — — - - - 
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Equations  for  the  higher  order  terms  in  (10)  and  (11),  can  also  be  obtained  but 
are  not  presented  here.  The  system  of  equations  becomes  closed  for  any  given 
cut-off  limit. 

s.  Solution  by  perturbation  expansion 

Instead  of  treating  the  system  of  highly  coupled  nonlinear  equations  (15)  “ 
(18)  directly,  we  shall  solve  them  by  expanding  the  dependent  variables  in  powers 
cf  a  coupling  parameter  a  so  as  to  separate  the  most  fundamental  part  F  (x  )  of 

oo 

Fq(x)  In  (17)  from  the  other  variables  to  make  F^  satisfy  an  uncoupled  equation, 
and  obtain  the  coupling  parts  of  the  various  functions  by  successive  approximation. 
The  expansions  which  serve  this  purpose  can  be  written  as 

Fo  "  Foo  +  aF01  +  0,2 F02  +  »  09) 


mo  -  raoo  +  <*"01  +  a  m02  +  ... 
F1  -aF10  +  a2FM  +a,Fi2  + 


iDj  ■  omjQ  +  a2mj  j  +  a*nij2... 
F2  "  0(2 F20  +  aSp21  + 


m2  ■  a2«n20  +  asro21  +  ... 


,  etc. 


Substituting  these  expansions  in  Eqs.  (15)  -  (18)  and  equating  the  coefficients 
of  the  various  powers  of  a  to  zero  we  then  obtain  the  following  system  of  equations; 
(first  set) 


(xF'*' 

oo 


+  (1+  F  )  F"  -  F'2]' 
v  oo7  oo  ooJ 


(20) 


IxF'j'q  +(1 


*  Foo>Fii  -  WioFio  +  5FiiF.ol' 


(20b) 
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(second  set) 


xm"  +  F 

m* 

oo  oo  oo 

xmji  +  F  mi¬ 
ll)  oo  10 

[xF 1 1  1  + 

L  01 

(1  + 

1 

m  m, 

2x2 

oo 

[xF'ii  ♦ 

(1  + 

0  , 

2F1  m,n 
oo  10 


.  -in 

10  oo 


(20c) 

(20d) 


oo 


oo 


10 


(21a) 


[xF'ii  +  <’  *  Foo>Fii  -  WioFi<  *  ^FooF1 1]  '  -  -  <F01Fii  -  *ilFio  +  3F'ilV 


OO 


+  <raio ♦  2moo"20>  ■ 6x2  ( > '  ■ 20F~  <  ^  > 


+  60^20  .  4of,20  _ 


xm‘ !  +  F  m' 
01  oo  01 


"  F01moo  "  2m10  » 


(21b) 

(21c) 


xra.'  i  +  F  »' .  -  2F 1  ro, , 
11  oo  11  oo  11 


"F0lm*  10  +  2F01ro10  "  3^F10rodl  +  F11moo)*,2m20  ' 


(21d) 


(Set  3»  limited  to  FQ2  and  mQ2  for  brevity) 


Hi  ♦  (l  ♦  Foo>Fii  -  2FooF02  +  FooF02]  1  ■  ■-  ("oo"1!!  *  molm10> 

*  [-FoiFii  +  ’oil1  -  "lo  <  jr  )■  -  6F0.  <  jr  >■  • 120  -r  ’  <22»> 

J  OO  O I 


xm* '  +  F  mi.  =  -2m, ,  -  F_ ,mi,  -  F__m' 

02  oo  02  II  01  01  02  oo 


(22b) 


Eqs.  (20a)  and  (20c)  are  found  in  Kuo  (1966,  1969)  while  (21b)  and  (21d)  differ 
somewhat  from  that  given  by  Kuo  due  to  the  different  6  used.  Observe  that  Eq.  (20a)  is 
nonlinear  but  it  contains  Fqq  alone,  hence  it  can  be  solved  together  with  the  boundary 
conditions  to  yield  F  ,  while  all  the  other  equations  are  linear  and  contain  lower 
order  functions  as  coefficients  and  non-homogeneous  terms.  Notice  also  that  the 
inclusion  of  the  Fj  term  in  (11)  makes  the  tangential  and  the  radial  velocities  coupled 
through  equations  (20d)  to  (22b). 
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i  -  0,1 
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o 

— 

i  =  0,1 

c.  The  Boundary  Conditions 

At  the  axis  we  require  that  each  of  the  variables  be  zero  for  the 
axisymmetric  flow,  i.e,,  at  x  =  0 

m  . 
oi 

m  a  • 

1 1 

At  the  rim,  x  “  xs>  we  assume  the  vertical  velocity  is  zero,  that  the 
radial  inflow  will  vary  slightly  with  height  while  the  rotation  rate  does 
not  vary  in  the  vertical.  Also,  we  shall  take  ail  the  radial  inflow  and 
rotation  in  the  first  terra  in  the  a  expansion,  and  none  in  the  higher 
perturbations.  Therefore  at  xg  the  solutions  satisfy 
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while  F  ,  Flrt  and  m  are  assigned  non-zero  values.  The  condition  on  the 
oo  lu  oo 

third  derivative  Is  needed  to  complete  the  requirement  of  four  boundary 
conditions  for  the  fourth  order  equations.  Integrating  (20a)  from  0  to  xg, 
applying  the  boundary  conditions  given  above  and  solving  for  F1 '  (x  )  yields 


oo 


F"  (x  )  -  F,io«»  -  F00<°> 
oo  s  . . — — 


(23) 


(24) 


(,+Foo(xs)) 

This  relation  shows  that  F"  (x  )  is  non-zero  if  F V*  (0)  i*  FjJ*  (0) .  A  similar 

OO  5  OO  OO 

procedure  applied  to  Eq.  (20b)  reveals  that 


(25) 
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F>0  <*5> 


Fio  <°> 


*FL.  (0)Fj0  (0) 


oo 


3F<^  <*s>F10<*s> 


(26) 


0 


+  F«>s» 

oo  s 


which  is  also  non-zero  if  the  numerator  Is  non-zero. 


At  the  rim  the  swirl  angle  6  can  be  written  in  terms  of  m  and  F 

o  o 

since  Eqs.  (12),  (14)  and  the  boundary  conditions  (24)  give  at  the  z*  «=  0  level 
0*-  tan  *(v/u)  =*  tan  '  (m  /2F  )  . 

o  O' 

Note  also  at  the.  rim  we  have  m  (x  )  *  m  (x  )  and  F  (x  )  **  F  (x  ). 

os  oo  s  os  oo  s 

d.  Details  of  the  Numerical  Procedure 

The  non-linear  ordinary  differential  equations  described  by  Eqs.  (20), 

(21)  and  (22)  with  their  corresponding  boundary  conditions  (23)  and  (24) 

are  solved  numerically,  in  the  order  listed  above,  using  a  shooting 

technique  (Conte  and  de  Boor,  1972).  Each  equation  is  rewritten  as  a 

system  of  first  order  equations  and  then  discretized  using  the  midpoint 

rule  (Kreiss  and  Oliger,  1973).  Thus,  e.g.,  each  fourth  order  equation  is 

reduced  to  a  system  of  four  first  order  equations  that  must  be  solved 

-4 

iteratively.  A  fine  grid  is  employed  near  the  axis,  e.g.,  Ax  *  10  , 

but  this  spacing  is  allowed  gradually  to  expand  to  a  much  coarser  net  near 
the  rim,  say  Ax  »  10  As  an  example,  a  total  of  575  grid  points  are  used 
on  the  interval  x  »  0  to  20.  Tests  using  nearly  double  this  number  of  grid 
points  yields  no  significant  differences. 

Utilizing  a  shooting  technique  to  solve  the  fourth  order  differential 
equations  for  the  various  F's  requires  that  Initially  a  guess  be  made  for 
the  first  and  second  derivatives  at  x  ■  0.  Using  these  guesses  the  third 
derivative  at  x  *  0  is  then  computed  directly  from  the  equation.  The 
system  of  four  first  c^der  equat'ons  are  solved  for  all  x  using  standard 
procedures  and  then  the  «•  «ue^  for  F}  an<j  jts  first  derivative  at  x  ■  x& 
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are  compared  with  the  desired  boundary  conditions.  If  the  differences  exceed 
a  specified  small  value,  say  10  **,  the  process  is  again  repeated  using  new 
estimates  to  replace  the  values  of  the  first  and  second  derivatives  at 
x  *  0.  These  new  values  are  obtained  using  Muller's  method  (Conte  and 
de  Boor,  1972). 

We  have  only  two  degrees  of  freedom  so  the  condition  on  the  third 

derivative  at  the  rim  can  not  be  strictly  enforced  but  we  except  only  those 

solutions  for  which  the  required  condition  is  naturally  approximately 

satisfied,  i.e.,  0,  i  =0,1,...,  J  =0,1,...  The  error  introduced 

*  »J 

by  this  procedure  is  quite  small  since,  e.g. ,  the  calculated  value  of  FJJ!  (x  ) 
varies  from  that  predicted  by  Eq.  (25)  typically  only  in  the  fourth  significant 
decimal  place. 

Our  equations  actually  possess  multiple  solutions  that  satisfy  all 

the  above  stated  boundary  conditions,  e.g.,  two  separate  solutions  for 

F  ,  Eq.  (20a),  are  known  to  exist  and  more  may  be  possible.  Using  initial 
oo 

guesses  with  F1  (0)  >  0  and  F* 1  (0)<  0  we  obtain  a  solution  which  contains 
oo  oo 

all  positive  values  of  F*  (one  cell  vortex)  while  the  second  solution, 

oo 

obtained  using  initial  guesses  F*  (0)<  0  and  F' 1 (0)  >  0,  contains  negative 

oo  oo 

values  of  F^  near  the  axis  and  positive  values  further  away  (two  cell 
vortices).  Kuo  (1967),  examining  a  vortex  in  an  unstably  stratified  atmos¬ 
phere,  shows  similar  two  cell  solutions  for  the  special  case  when  all  the 
derivatives  are  zero  at  the  rim.  Our  two  cell  solutions  remain  essentially 

i 

unchanged  near- the  axis  as  the  rim  is  moved  further  and  further  away.  This 

behavior  is  not  observed  in  the  laboratory  vortex  (Ward,  1972  and  Church 

et  al.,  1979),  so  in  this  study  we  disregard  this  solution  and  utilize  only 

initial  guesses  that  satisfy  F'  (0)  >  0  and  F1'  (0)  <  0.  The  exact  magnitude 

oo  oo 

for  these  choices  depends  on  the  strength  of  the  inflow  at  the  rim. 

The  solution  of  the  m  equation  is  straight  forward  and  does  not  change 
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significantly  with  the  change  of  the  vertical  velocity  at  the  axis  and  there¬ 
fore  it  will  not  be  presented  here. 

4.  Discussion  of  Results 

In  this  study  a  similarity  approach  is  utilized  to  obtain  the  solution  of 
the  two-dimensional  nonlinear  vortex  equations  for  comparison  with  the  results 
obtained  from  laboratory  vortex  simulators,  so  that  the  vertical  and  radial 
variations  are  uncoupled.  For  simplicity,  the  geometry  of  the  vortex  simulator 
is  ignored  and  only  the  nature  of  the  flow  within  the  convergence  or  inflow 
region  is  investigated  for  a  range  of  swirl  angles.  Care  must  be  exercised 
in  interpreting  our  results  however  since  the  similarity  transform  ignores 
the  geometry  and  by  limiting  this  study  to  within  the  convergence  region  we 
also  ignore  the  source  of  the  convergence,  i.e.,  the  exhaust  fan,  and  thus 
there  is  no  guarantee  that  different  rim  rotation  rates,  with  the  same 
specified  radial  inflow,  will  experience  the  same  volume  flow  rate.  Fortunately 
Davies- Jones  (1973)  has  shown  that  within  the  experimental  ranges  for  the 
three  nondimens ional  quantities,  namely,  the  radial  Reynolds  number,  the 
aspect  ratio  and  the  swirl  ratio, it  is  the  swirl  ratio  that  most 
strongly  controls  the  laboratory  vortex  flow  configuration.  In  our  calculations 
some  small  variations  of  the  aspect  ratio  and  radial  Reynolds  number  are  inherent 
because  of  the  nature  of  the  similarity  approach  but  these  changes  are  thought 
to  be  minimal  and  should  not  distort  the  general  conclusions.  Also  our  model 
includes  the  viscous  or  diffusion  terms  but  does  not  include  the  boundary  layer 
explicitly  since  the  motion  predicted  at  z*  =  0  is  non-zero  and  not  the  no-slip 
flow  required  in  boundary  layer  theory.  The  zero  value  for  z*  should  be 
interpreted  as  occurring  near  the  top  of  the  surface  boundary  layer. 

We  will  vary  the  boundary  conditions  at  the  rim  (x  ■  xg)  representing 
changes  in  the  strength  of  the  radial  velocity  or  inflow  (Fqo(xs)),  the 
vertical  variation  of  the  radial  velocity  component  (fiq(xs) *F20^xs^  ’  " 


and  the  strength  of  the  rotation  or  tangential  velocity  (m00(xs))*  and 
observe  the  predicted  flow  pattern.  Note  that  the  height  dependency  is 
determined  by  the  first  and  higher  order  terms  in  Eq.  (II)  and  is  established 
by  the  choice  of  boundary  conditions  at  the  rim.  Results  will  be  presented 
normally  for  two  boundary  parameters  fixed  while  the  third  is  varied  over 
some  range.  The  resulting  flows  may  not  always  be  observed  since  in  the 
laboratory  vortex  there  may  be  some  natural  compensation  at  the  rim  in  the 
vertical  variation,  for  example,  as  changes  occur  there  in  either  the 
strength  of  the  radial  inflow  or  rotation  rate.  Nevertheless  we  can  clearly 
show  the  trends  produced  by  individual  changes  in  these  boundary  parameters. 

Our  results  illustrated  in  Figs.  I  through  11  are  displayed  in  terms  of 

1  /2 

the  dimensional  radial  variable  r  -  r*  or  in  (r*)  scale  to  better  visualize 

the  predicted  flow  pattern  near  the  axis.  The  first  10  figures  are  for  an 
intense  vortex  similar  to  those  produced  in  modern  laboratory  vortex 
simulators.  The  outer  rim  is  located  at  r*(x&)  *  1.932  m  while  the  last 
figure  shows  a  weaker  configuration  with  maximum  radius  of  0.642  m.  Be¬ 
cause  the  intensity  of  the  inflow  and  rotation  at  the  outer  rim  greatly 
effects  the  location  of  the  maximum  tangential  velocity  in  the  x  coordinate 


an  appropriate  choice  for  rQ  must  be  established  for  each  of  the  two  cases 
to  convert  our  calculations  in  x  back  Into  r  ■  r*  scale.  The  scaling  factor 

r  2  t  /^r  j “2 

r  is  calculated  once  for  each  case  via  the  formula  x. _  *=  vmax  *  '  o' 

o  vmax 

by  selecting  0.15  m  as  the  appropriate  radius  for  the  maximum  tangential 

velocity,  in  approximate  agreement  with  laboratory  simulations  (Church 

et  al . ,  1979).  Both  radii  given  above  correspond  to  a  maximum  value  for 

x  of  20.  This  value  is  selected  as  typical  from  a  large  number  of  tests. 

Using  numerical  results  representative  of  each  case  the  scaling  factor  for 

the  more  intense  flow  is  determined  as  r  =  0.215  m  while  the  weaker  case  uses 

o 

r^  ■  0.0728  m.  The  equations  for  Fq,  mQ,  ...,  themselves  are  independent 
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of  rQ,  6  and  v  but  to  convert  back  to  the  dimensional  variables  requires  some 
appropriate  choices.  The  radial  distance  is  labeled  r*  in  the  figures  indicating 
that  our  choice  of  rQ  has  been  used  to  convert  from  the  x  coordinate  back  into 
the  dimensional  r  =  r*  coordinate. 

Radial  profiles  of  P  and  its  components  are  shown  in  Fig.  1.  These  are 
computed  using  F  (x  )  *  250,  F._(x_)  =  -  0.8  and  m  (x  )  =  19*»  and  are  presented 

OO  5  I U  5  OO  S 
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in  a  uniform  (r*)  scale.  F^  and  its  components  should  be  Interpreted  as  a 

measure  of  the  vertical  velocity  (Eq.  13).  The  solid  curve  labeled  F*  ,  obtained 

oo 

from  (20a)  has  its  maximum  value  at  the  axis,  i.e.  r*  =  0,  and  by  itself  is  typical 
of  the  profile  found  in  a  one-cell  vortex  (Kuo,  1966).  When  higher  order  expansion 
terms,  obtained  from  Eqs.  (21a)  and  (22a)  are  included  in  F^  *  F^  +  aF^j  +  a2F^2 
with  a  taken  as  one,  the  profile  is  radically  different  since  F^  =  0  at  r*  =  0 
and  thus  represents  the  initial  stage  of  a  two  cell  vortex.  The  latter  is  generally 
characterized  by  negative  vertical  velocities  at  and  near  the  axis  with  non-negative 
values  further  away.  For  our  purposes  three  terms  will  be  Sufficient  to  approximate 
Fq  or  F^,  since  e.g.,  the  dotted  curve  representing  Fg.,  is  significantly  smaller  in 
magnitude  than  either  F^q  or  F^  indicating  the  decreasing  importance  of  higher 
order  terms  in  the  a  expansion.  Higher  order  terms  in  the  (6z*2)power  series 
expansion  Eqs.  (10)  and  (11)  like  Fj ,  where  Fj  =  aF^g  +  a2Fjj,  are  similarly  much 
smaller  than  the  zeroth-order  terms.  Note  that  lO.Fjg  as  plotted  in  Fig.  I  is 
still  small  In  magnitude.  The  second  term  in  Fj ,  i.e.  F^ j ,  is  even  smaller.  We 
will  neglect  the  very  small  contributions  of  F£  and  other  higher  order  terms. 

Included  in  Fig.  1  is  the  dashed-dot  curve  showing  the  radial  profile  of  F^ 
when  terms  generated  by  vertical  diffusion  are  removed  from  the  calculations.  The 
resulting  differences  in  F^  are  clearly  small  as  seen  by  comparing  the  dashed-dot 
and  solid  curves  and  do  not  change  the  two  cell  nature  of  the  vortex  which  results 
when  radial  and  tangential  motions  are  coupled.  Our  calculations  show  that  as  the 
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vortex  becomes  more  intense  this  small  contribution  is  reduced  even  further  (not 

displayed)  while  in  the  convergence  region  of  a  weak  vortex  vertical  variations 

are  more  pronounced  so  that  vertical  diffusion  becomes  more  important,  as  shown 

in  Fig.  11.  In  this  weaker  vortex  configuration,  Fig.  11,  the  ratio  of  F,_/F 

i  u  oo 

at  xs  is  an  order  of  magnitude  larger  than  that  used  to  compute  the  more  intense 
vortex  described  in  Figs.  1  through  10. 

According  to  Hail  (1972)  the  vortex  core  is  characterized,  above  the  boundary 
layer,  by  a  slight  spreading  out  of  the  core  with  height  hence  a  small  decrease 
in  radial  and  tengential  velocities  with  height  and  the  development  of  an  adverse 
pressure  gradient. 

We  tune  our  model  to  produce  similar  behavior  by  choosing  Fj  as  minutely 

negative  at  the  rim,  hence  negative  over  the  entire  radius,  giving  an  increasingly 

negative  contribution  with  increasing  height  to  the  radial  velocity  component  and 

consequently  to  the  tangential  velocity.  Fig.  2  shows  that  the  value  of  F^  is 

sensitive  to  changes  in  Fj  since  non-zero  Fj  values  allow  the  tangential  velocity 

and  radial  motions  to  be  coupled.  If  F,  is  identically  zero  everywhere  then  F 

i  o 

wiil  be  identically  F^,  i.e.,  the  radial  motion  becomes  independent  of  the  rota¬ 
tion  rate.  The  magnitude  of  F^  is  chosen  so  that  the  coupling  effects  of  Fj  and 

ml  produce  nearly  zero  values  at  the  axis  for  F^,  solid  line,  when  the  swirl 

angle,  0  *  tan  *  m  /2F  ,  is  in  the  neighborhood  of  twenty  degrees  in  approximate 
o  o 

agreement  with  results  found  in  laboratory  simulations  (Church  et.  al.,  1979, 
Church  and  Snow  1979). 

In  Fig.  3,  radial  profiles  of  F^  are  presented  for  four  different  rotation 
rates  at  the  rim,  I.e.,  with  moQ(xJ)  ■  170,  1 94 ,  2^0  and  280  while  F^fx^  ■  250 
and  Fjg(xg)  “  “  0.8  are  each  held  fixed.  This  is  equivalent  to  varying  the  swirl 
angle  at  the  rim  from  18.78°  to  29.25°  (z*  -  0)  by  adjusting  the  tangential  velo- 

i 

t 

j  city.  It  should  be  noted  that  laboratory  simulations  measure  their  swirl  angle 

t 
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somewhere  between  the  rotating  screen  and  the  updraft  core  and  not  at  the  rim 

exactly.  For  the  smaller  value  of  m00(xs)  of  170  the  F^  profile,  given  by  the 

dashed  curve,  is  positive  throughout  the  entire  radius  but  already  significantly 

reduced  near  the  axis  as  seen  when  compared  with  the  solid  curve  labeled  F1 

oo 

which  does'not  include  the  coupling  effect  of  the  tangential  flow.  As  the  rota¬ 
tion  rate  is  increased  F^  becomes  increasingly  smaller  at  the  axis  and  is  nearly 
zero  when  nioo(xs)  =  194  (solid  line).  Further  increases  in  the  rotation  rate 
produce  negative  F^  values  at  the  axis  typical  of  the  two  celled  vortex.  Note 
that  the  profile  away  from  the  axis,  say  r*  >  0.25  m,  remains  almost  uneffected 
by  variations  in  the  swirl  angle.  In  the  laboratory  simulator  where  the  volume 
flow  rate  is  unchanged  in  the  experiment  there  would  be  somewhat  more  of  an 
increase  in  this  outer  region,  of  the  order  of  10%,  to  compensate  for  the  down- 
watf  motion  at  the  axis.  Thus  our  results  for  increasing  values  of  mQ  at  the 
boundary  have  a  small  decreasing  volume  f low  rate  for  the  same  radial  inflow 
rate.  This  may  be  interpreted  as  having  a  slightly  decreasing  aspect  ratio  and 
radial  Reynolds  number.  Consequently  the  swirl  ratio  S  will  increase  more  than 
that  predicted  due  to  changes  in  the  swirl  angle  alone  since  S  is  inversely  pro¬ 
portional  to  the  aspect  ratio. 

-1/2 

In  Fig.  4  values  of  x  F  ,  which  is  proportional  to  the  radial  velocity, 

are  given  up  to  a  radius  of  0.28  m.  The  two  cell  flow  is  very  evident  in  the 

dashed-dot  (m  =  240)  and  dotted  curves  (m  „  ■  280)  since  both  have  negative 
oo  oo 
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values  next  to  the  axis.  The  solid  curve  labeled  x  Fjq  shows  the  smallness 
of  the  negative  first  order  term  when  compared  with  zero  order  terms. 

The  development  of  the  downward  vertical  motion  at  the  axis  is  clearly  seen 
to  descend  from  high  levels  of  z,  in  agreement  with  laboratory  findings  (Church 
et  al .  1979)*  since  negative  vertical  motion  is  possible  when  small  positive 
zeroth-order  and  negative  first-order  terms  are  combined  provided  at  some  z*  >  0 


the  expression  |6z*2Fj|  >  P  in  Eq.  (13)  is  first  satisfied.  This  occurs  when 
F^  is  much  subdued  but  still  positive  at  the  axis,  e.g.,  as  does  occur  as  shown 
in  Fig.  3*  The  combination  of  zeroth  and  first  order  terms  would  then  lead  to 
a  reduction  with  height  in  the  vertical  velocity,  a  stagnation  or  zero  point 
and  finally  negative  vertical  motion  above  that  level.  This  type  of  motion  is 
clearly  seen  in  the  laboratory  simulations  of  the  vortex  (Church  et  al . ,  1979, 

Ward  1972). 

Profiles  of  the  tangential  velocities  are  shown  in  Fig.  5  for  the  same  cases 

discussed  in  Figs.  3  and  4.  As  the  swirl  angle  at  z*  =  0  is  increased  in  value 

for  four  cases  between  18.78  to  29.25°,  corresponding  to  the  dashed  through 

dotted  curves  respectively,  the  tangential  velocity  maximum  moves  further  from 

the  axis,  a  feature  commonly  observed  in  laboratory  simulations  (Church  et  al. , 

1979).  Outside  the  radius  of  maximum  tangential  flow  the  shapes  of  the  profiles 

are  very  similar.  This  follows  since  the  value  of  the  zeroth  order  angular 

momentum,  mQ  ■  vQr,  remains  nearly  constant.  For  example  with  a  value  of  194  at 

the  rim  the  value  of  mQ  is  192.8  at  r*  *  0.2422m  which  is  very  close  to  the  radius 

of  maximum  velocity  at  0.15m.  This  means  that  the  vertical  velocity  is  confined 

within  a  region  certainly  less  than  twice  the  radius  of  the  maximum  tangential 

velocity.  The  flow  behavior  very  near  the  axis  is  in  a  state  of  near  solid 

rotation  thus  Fq,  mQ,  ...  each  are  nearly  proportional  to  x  as  shown  in  Table  l 

for  the  case  m  (x  )  ■  194.  Obviously  the  slope  there  will  vary  as  the  swirl 
OO  s 

angle  is  changed.  Note  that  as  mQ(xs)  is  increased  from  170  to  194,  dashed  and 
solid  line  respectively  in  Fig.  5,  the  slppe  of  the  tangential  velocity  is 
slightly  reduced  near  the  axis  but  enhanced  as  it  approaches  its  maximum  value. 
Further  increases  In  the  swirl  angle  continue  to  steepen  the  slope  just  inside 
the  radius  of  maximum  velocity. 

The  tangential  velocity  is  observed  to  become  negative  in  the  region  of  the 
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core  where  the  two  cell  vortex  contains  negative  vertical  velocity  and  outward 
radial  velocity  and  hence  gives  rise  to  an  advective  loss  of  angular  momentum 
in  the  core.  In  the  laboratory  simulator  the  downward  motion  at  the  axis  in 
the  two  cell  vortex  appears  to  descend  downward  through  the  baffling  near  the 
top  of  the  vortex  generator  down  to  the  lower  surface.  Obviously  our  model 
which  is  restricted  to  the  convergence  region  of  the  vortex  cannot  duplicate 
this  behavior.  The  configuration  given  by  the  dashed-dot  and  dotted  curves  in 
Fig.  5,  representing  predicted  rotational  flow  for  a  two  celled  vortex,  will 
give  rise  to  inertial  instability  since  the  gradient  of  the  circulation  changes 
signs  (Rayleigh,  1916  and  Synge,  1938). 

Another  type  of  instability  is  revealed  in  Fig.  6  by  the  radial  profiles 

of  m^,  here  nT  is  proportional  to  the  vertical  vorticity  and  the  prime  denotes 

a  derivative  with  respect  to  x.  Even  for  the  smaller  boundary  value  of  m  (xg)  = 

170  (dashed  curve)  the  gradient  of  m^  changes  sign  revealing  that  the  flow 

configuration  is  barotropical ly  unstable.  The  slope  of  the  uncoupled  or  zeroth 

order  term  m*  is  of  the  same  sign  thus  the  instability  portrayed  in  other  curves 
oo 

arises  because  of  the  coupling  effect  of  the  radial  and  vertical  flows  back  on 
the  tangential  motion.  Note  that  the  maximum  value  for  m^  occurs  at  the  axis 
while  in  the  other  curves,  representing  for  increasing  swirl  angles,  the 
maxima  are  achieved  progressively  further  away  from  the  axis  and  closer  to  the 
tangential  velocity  maxima.  In  fact  these  maxima  occur  where  the  gradient  of 
the  tangential  velocity  is  very  large  as  seen  when  comparing  Figs.  5  and  6. 

These  results  indicate  that  azimuthal ly  varying  three-dimensional  disturbances 
precluded  from  the  present  axisymmetric  model  will  be  created  under  these  unstabl 
conditions,  which  will  even  out  the  negative  vorticity  and  negative  rotation  in 
the  core  to  make  them  not  observable  in  laboratory  simulation.  The  possibility 
that  the  barotropic  instability  of  the  tangential  flow  profile  may  lead  to 


generation  of  suction  vortex  type  disturbance  in  tornadoes  as  suggested  by 
Fujita  (1972)  has  been  investigated  by  Staley  and  Gall  (1979).  Snow  (1978)  has 
postulated  the  other  alternative,  i.e.,  an  inertial  instability. 

Table  I.  Values  of  F  ,  F.,  m  ,  m.  and  their  components  at  two  locations  near 

o  1  o  I 

the  axis.  All  F  and  m  values  are  to  be  multiplied  by  10 


0.0001  0.0043  196  -146 

0.001  0.01359  1976  -1453 


-49 

-499 


-0.9  13736 

-8.6  109742 
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The  radial  profile  of  the  scaled  velocity  magnitude  (|v|v  .10  )  at 

level  2  -  0  and  its  tangential  and  radial  component  are  represented  by  the 
continuous,  the  dashed  and  the  dotted  curves  in  Fig.  7  respectively.  This 
form  of  presentation  is  taken  since  a  conversion  to  actual  velocity  would 
require  the  knowledge  of  the  viscosity  coefficient  v.  Our  model  predicts 
less  than  a  four-fold  increase  in  the  velocity  magnitude  from  the  rim  to  its 
maximum  value.  This  appears  to  agree  very  well  with  laboratory  simulations 
(Church  etal . ,  1979).  Also  note  that  the  radial  component's  contribution 
to  the  velocity  magnitude  is  very  small  inside  the  radius  of  maximum  tengen- 
tial  velocity  but  dominates  near  the  rim. 

For  very  small  swirl  angles  (not  shown),  our  model,  in  agreement  with 
laboratory  findings  (Snow  et  al . ,  1980)  gives  weakly  swirling  flows  without 
a  central  core.  In  this  case  the  tangential  component  increases  with  radius 
right  up  to  the  rim. 

Radial  profiles  proportional  to  the  vertical  velocity,  the  tangential 
velocity  and  the  vertical  vorticity  are  shown  in  Fig.  8  through  10,  respective¬ 
ly.  Three  different  values  of  the  swirl  angle  ranging  from  17-92  to  23-32° 

(z*  =0)  are  used  by  varying  the  radial  inflow  component  while  keeping  the 

rotation  rate  at  the  rim  constant.  As  the  boundary  value  of  F  is  decreased 

oo 

from  300  (dashed  curve)  to  225  (dotted  curve)  while  m^Xg)  *  19^  and 
Fi0(xs)  *  -  0.8  are  each  held  constant,  the  vertical  flow  is  found  to  change 
from  all  positive  values  to  the  two  cell  configuration,  in  agreement  with  what 
happens  as  the  swirl  angle  is  increased.  Correspondingly,  the  tangential 
velocity  (Fig.  9)  shows  a  decrease  in  magnitude  with  an  outward  expansion  of 
the  radius  of  the  maximum  value  in  accordance  with  the  law  of  conservation  of 
angular  momentum.  The  curves  in  Fig.  10  show  a  similar  pattern  to  that 
discussed  earlier  for  the  vertical  vorticity. 


Our  equations  are  non-linear  so  the  nature  of  their  solutions  might  vary 

somewhat  as  the  magnitudes  of  the  boundary  values  are  varied  but  we  find  that 

fora  large  range  of  boundary  conditions  they  still  behave  similarly.  A 

weaker  flow  configuration  using  F  (x  )  a  25  is  shown  in  Fig.  11,  but  even 

oo  s 

in  this  case  it  still  is  possible  to  have  the  vertical  velocity  change  from 
all  positive  (dashed  line,  m00(*s)  “  10)  to  that  found  in  a  two-celled  vortex 
(dotted  line,  mOQ(xs)  =  30)  by  changing  the  tangential  velocity  at  the  rim. 

No  negative  tangential  velocities  are  generated  for  this  case  where  the  swirl 
angle  is  varied  from  11. 31  to  30  .96°.  Tests  using  F^  (x&)  »  500  also 
generate  the  two-celled  vortex  with  the  appropriate  choices  of  m00(xs)*  Note 
that  when  our  results  are  converted  into  actual  velocities  the  choice  for  the 
value  of  viscosity,  used  in  scaling  the  equation,  does  not  determine  the  nature 
of  the  flow  since  the  same  flow  characteristics  are  found  over  a  broad  range 
of  boundary  values. 

5.  Conclusion 

In  this  study  a  model  of  the  convergence  region  of  a  laboratory  vortex 
is  investigated  via  a  similarity  approach.  The  predicted  flow  is  shown  to 
vary  significantly  only  near  the  axis  as  the  swirl  angle  at  the  outer  rim 
is  changed.  This  type  of  response  appears  to  be  exactly  similar  to  that 
observed  in  the  laboratory.  Of  course,  not  all  the  details  of  the  flow 
found  in  the  laboratory  vortex  generator  can  be  reproduced  here,  but  the 
nature  or  trend  of  changes  predicted  in  this  study  appears  to  be  generally 
val id. 
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List  of  Figures 

Fig.  1.  Radial  profiles  of  F^  which  is  proportional  to  the  vertical 

velocity,  and  its  components  are  given.  Fjq  is  also  displayed. 

First  subscript  denotes  the  6 z2  expansion  while  the  second  Indicates 

location  in  a  expansion.  Boundary  values:  F  =  250,  m  **  194 

oo  oo 

and  Fjq  =  -0.8.  Dashed-dot  curve  is  F^  when  vertical  diffusion  is  ignored. 

Fig.  2.  Radial  profiles  of  F^  computed  using  three  different  boundary 

values  for  F^q,  i.e,  (solid)  Fjq  =  -0.8,  (dotted)Fjg  =-l.l  and 

(dashed)  Fjg  -  -0.6.  Other  boundary  values  same  as  in  Fig.  1. 

Fig.  3.  Radial  profiles  of  F^  computed  for  four  different  swirl  angles 

obtained  by  varying  the  tangential  velocity.  Boundary  values  used  are 

(solid)  m  -  194 ,  (dashed)  rann  =  170,  (dashed-dot)  m  =  240  and 
OO  Uw  oo 

(dotted)  mQo  =  280.  Other  boundary  values  same  as  in  Fig.  1. 

- 1/2 

Fig.  4.  Same  as  Fig  3  except  showing  x  Fq,  which  is  proportional  to 

the  radial  velocity,  on  a  linear  scale  near  the  axis. 

-1/2 

Fig.  5.  Same  as  Fig.  3  except  showing  x  mQ  which  is  proportional 
to  the  tangential  velocity. 

Fig.  6.  Same  as  Fig.  3  except  showing  which  is  proportional  to  the 
vertical  vortlclty. 

-2  1/2  -1  “1 
Fig.  7-  Radial  profiles  (time  10  )  of  (u2+v2)  v  (solid),  w 

(dashed)  and  uv  ^(dotted).  Boundary  values  same  as  used  in  Fig.  1. 

Fig.  8.  Radial  profiles  of  F^  for  three  swirl  angles  obtained  by  varying 
the  radial  velocity  component.  Boundary  values:  (solid)  F  ■  250, 

OO 

(dashed)  F  ■  300  and  (dotted)  F _  ■  225.  Other  boundary  values 

OO  oo 

same  as  Fig.  1 . 
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Fig.  9.  Same  as  Fig.  8  except  showing  x  m  . 

o 

Fig.  10.  Same  as  Fig.  8  except  showing  m^. 

Fig.  11.  Radial  profiles  of  F^  for  three  swirl  angles.  Computed  using 

•  boundary  values  (solid)  m  =  20,  (dashed)  m  =10  and  (dotted) 

00  00 

m  =  30.  Otherwise  ^0Q(XS)  =  25  and  ^jq(xs)  =~0.3.  The  dashed-dot 
curve  is  F'  when  vertical  diffusion  is  ignored. 
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